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Preface 

It  is  interesting  to  note  that  only  a  short  decade  or  sc  ago, 
neither  the  theory  nor  the  computational  capability  (let  alonr  the 
requirement)  for  attacking  the  problem  of  this  thesis  topic  existed. 

In  the  early  1960's,  ?..  E.  Kalman  presented  his  theories  on  optimal 
control  and  filtering  technioues.  The  advent  of  that  electronic  marvel 
of  our  age,  the  digital  computer,  provided  the  tool  for  applying 
Kalman's  filtering  theory.  Many  Saturday  and  Sunday  mornings  as  veil 
as  veeknlghts  were  spent  with  our  dear  friend,  the  CDC  6600  computer, 
in  order  to  bring  to  you,  the  reader,  the  meticulous  CALCOMP  plots  of 
the  system  state  variable  covariances  presented  throughout  this  report. 

Although  there  were  many  individuals  from  both  the  Air  Force 
Avionics  Laboratory  and  the  School  of  Engineering  faculty  who  rendered 
us  assistance  from  time  to  time  during  the  preparation  of  this  thesis, 
there  are  two  people  to  whom  we  are  especially  indebted.  We  would  like 
to  express  our  most  sincere  appreciation  to  Lt  Peter  S.  Naybeck,  faculty 
advisor  for  this  thesis  topic.  Even  though  he  was  enduring  the  trials 
and  tribulations  of  an  expectant  father  and  also  attending  the  Academic 
Instructor  School  (in  addition  to  his  everyday  faculty  duties)  during 
the  time  period  of  this  study,  he  somehow  always  managed  to  find  time 
to  talk  with  us.  For  his  patient  guidance  and  instruction  we  are  indeed 
grateful.  We  also  gratefully  acknowledge  the  assistance  of  Mr.  Richard 
M.  Reeves  of  the  Avionics  Laboratory  V.  providing  us  with  computer 
programs  and  other  documents  related  n  our  study.  We  especially  . 
thank  him  for  sharing  with  us  hi*;  «•:  >  rcc  and  expertise  in  the 

application  of  Kalman  filter  design  tm-ory. 
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Finally,  as  we  emerge  from  Che  dark  "Caves  of  Knowledge"  blinking 
our  eyes  in  response  to  the  bright  sunlight,  we  would  like  to  express 
appreciation  to  our  wives,  Dottie  and  Carolyn,  who  arc  pAtiently  waiting 
to  welcome  us  hone. 


Ronald  R.  Butler 
Ceorge  T.  Rhue 
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s  Abstract 

This  report  Is  a  Kalman  filter  design  study  for  the  proposed 
Integrated  Navigation  itellite/Inertlal  System  (INI).  Primary 
emphasis  is  placed  upon  determination  of  the  "best"  filter  state 
variable  vector  and  investigation  of  various  measurement  rates  using 
external  range  measurements  from  a  set  of  27  non-synchronous  sate¬ 
llites.  The  INI  system  errors  are  assumed  to  be  represented  by  a 
44  state  linear  system  model,  and  the  filter  operates  without  benefit 
of  an  altimeter.  A  one-hour  INI  flight  at  constant  speed  and  altitude 
over  a  great  circle  path  is  simulated  on  the  digital  computer  and 
filter  designs  are  compared  by  plotting  the  system  position,  velocity, 
and  attitude  error  covariances  versus  tine.  A  15  state  filter  with 
weak  coupling  terms  removed  is  determined  to  provide  the  best  tradeoff 
between  accuracy  and  computational  burden.  Filter  performance  is 
compared  at  5,  15,  30,  60  and  90  second  measurement  update  rates. 
Additionally,  "optimal"  sequencing  of  satellite  observables  is  shown 
to  provide  improved  performance. 
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KALMAN  FILTER  DESIGN  FOR  AN  INERTIAL 
NAVIGATION  SYSTEM  AIDED  BY  NON-SYNCHRONOUS 
NAVIGATION  SATELLITE  CONSTELLATIONS 

I.  Introduction 

This  report  Is  a  Kalman  filter  design  study  for  the  proposed  Inte¬ 
grated  Navigation  Satellite/ Inertial  System  (INI).  This  design  is  a 
aubproblem  of  the  proposed  INI  system  which  Is  a  current  research  and 
development  topic  of  interest  to  the  Navigation  Division  of  the  Air 
Force  Avionics  Laboratory  at  Wrlght-Patterson  Air  Force  Base  as  well  as 
other  government  and  private  agencies. 

The  development  of  this  study  is  presented  by  chapters  in  the 
following  sequence.  The  preliminaries  such  as  reference  frames  and 
their  angular  rates  are  presented  in  the  second  chapter.  Chapter  III 
discusses  theory  which  is  peculiar  to  the  use  of  navigation  satellites 
as  an  updating  aid  to  inertial  navigation  systema.  Chapter  IV  is  a 
development  of  the  Kalman  filter  equations  required  to  solve  the  naviga¬ 
tion  problem.  The  forty-four  state  system  model  selected  for  this  study 
is  presented  next  in  Chapter  V.  The  intent  of  Chapter  VI  is  to  deter¬ 
mine  the  "best"  state  variable  vector,  both  In  terms  of  dimension  and 
selection  of  components  for  implementation  in  the  filter.  Chapter  VII 
investigates  the  effects  of  varying  the  measurement  update  rate  of  this 
design.  A  description  of  the  computer  program  used  in  simulating  the 
INI  flight  is  given  in  Chapter  VIII.  The  results  and  eoncluaions  com¬ 
pose  the  final  chapter  of  this  design  study. 

In  tearated  NAVSAT/Inertial  System 

The  basic  idea  of  the  INI  is  to  combine  external  navigation 
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information  (such  as  range,  range  rate,  and  attitude)  received  from 
orbiting  constellations  of  satellites  with  the  inertial  navigation 
system  (INS)  of  an  aircraft  or  missile  to  provide  highly  accurate  posi¬ 
tion  and  velocity  information,  such  as  would  be  required  for  Instrument 
landing  or  weapon  delivery.  The  Navigation  Satellite  (NAVSAT)  portion 
of  the  system  consists  of  a  total  of  27  orbiting  satellites.  The  num¬ 
ber  of  satellites  and  their  orbits  are  arbitrary  and  the  particular 
configuration  selected  for  this  study  will  be  discussed  in  detail  in 
Chapter  III.  A  filter  design  study  utilising  a  cluster  of  four  satel¬ 
lites  in  a  Y-conflguratlon  with  the  central  satellite  in  synchronous 
orbit  and  the  remaining  satellites  orbiting  about  it  is  detailed  xn 
Reference  1  of  the  bibliography.  There  are  several  existing  satellites 
under  consideration  for  adaptation,  one  of  which,  the  Tlmation  III,  will 
be  used  as  a  baseline.  Each  satellite  contains  a  transmitter,  receiver, 
and  a  clock  (i.e.,  quarts  crystal  oscillator).  A  ground  tracking  net¬ 
work  periodically  measures  and  updates  the  ephem«.ris  and  satellite  clock 
phase  and  frequency  so  as  to  maintain  synchronisation  of  all  satellite 
clocks.  The  satellites  in  turn  continually  transmit  this  information 
together  with  identifiable  range  codes  to  the  user.  The  signals  from 
each  satellite  are  modulated  by  orthogonal  codes  so  that  they  may  be 
distinguished  from  each  other  by  the  user.  By  means  of  a  correlation 
detector  the  time  shift  between  each  satellite  signal  and  the  user's  un¬ 
synchronized  clock  is  measured  in  the  user's  receiver  (Ref  1:1-2  &  1-3). 
The  portion  of  the  system  on  board  the  aircraft  consists  of  a  computer, 
receiver,  and  an  inertial  measurement  unit  (IMU).  The  receiver  is  re¬ 
quired  for  acquisition  of  the  satellite  signal.  The  IMU  consists  basi¬ 
cally  of  a  set  of  accelerometers  and  gyros  which  provide  internal  or 
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on-board  navigation  information  using  specific  force  and  angular  rates 
along  three  mutually  orthogonal  axes.  The  computer  is  used  to  apply  the 
Kalman  filter  equations,  i.e,,  to  compute  the  Kalman  gains  or  weighting 
coefficients,  to  estimate  the  state  variables,  and  to  apply  the  control 
to  the  Ih'S. 

Kalman  Filter 

The  Kalman  filter  is  simply  an  optimal  recursive  data  processing 
algorithm  which  resides  in  the  central  processor  or  on-board  computer. 
This  filter,  or  computer  program,  combines  all  available  measurement 
data,  plus  prior  knowledge  of  the  system  and  measuring  devices  to  pro¬ 
duce  an  estimate  of  the  desired  variables  in  such  a  manner  that  the 
resulting  error  Is  minimised  statistically*  Stated  in  simpler  terms, 
it  provides  the  best  estimate  possible  subject  to  certain  modeling 
assumptions. 

The  filter  will  enhance  the  attitude  and  position  information 
accuracy  by  weighting  each  data  source  heavily  in  the  frequency  regime 
where  it  provides  good  information,  and  suppressing  it  in  the  region 
where  It  is  in  error.  The  inertial  system  provides  good  high  frequency 
Information,  but  it  drifts  slowly  and  therefore  exhibits  poor  low  fre¬ 
quency  performance.  On  the  other  hand,  the  NAVSAT  data  is  good  on  the 
average,  but  subject  to  high  frequency  noise.  Thus  the  filter  will  use 
the  good  low  frequency  NAVSAT  information  to  damp  out  the  slowly  growing 
errors  Inherent  in  the  inertial  systc.  . 

Basic  Assumptions.  Under  the  t'a.-.u  restrictions  of  system  linear¬ 
ity,  noise  whlteneos,  and  Gaussirmc  ••  noises,  the  Kalman  filter  can 
be  shown  to  be  the  best  filter  of  ...v  . ...ccivuble  form.  Although  the 


3 


GA/EE/74M-1 


system  Itself  is  actually  nonlinear ,  the  formulation  of  an  approximate 
linear  error  model  makes  linear  analysis  possible.  The  justifications 
for  the  linear  system  model  are  that  (1)  the  use  of  linear  models  in 
engineering  studies  has  proved  fruitful  and  (2)  the  techniques  of  linear 
systems  analysis  are  well  developed,  whereas  those  for  nonlinear  systems 
are  not,  in  general.  "Whiteness"  Implies  that  the  noise  value  Is  not 
correlated  In  time  and  also  has  equal  power  at  all  frequencies. 
Gauesienness  pertains  to  the  amplitude  of  the  noise;  at  any  single  point 
In  time,  the  probability  density  of  a  Gaussian  noise  amplitude  takes 
on  the  shape  of  a  normal  bell-shaped  curve.  These  three  assumptions 
greatly  simplify  the  mathematics  of  the  problem  and,  In  fact,  render 
them  tractable  (Ref  2:3-11  &  Ref  3:132). 

Type  of  Filter  Employed  In  the  INI.  An  Indirect  filter  estimates 
the  errors  in  the  navigation  information  using  the  difference  between 
IRS  and  external  NAVSAT  data  as  the  measurements  to  drive  the  filter. 

This  is  In  contrast  to  the  direct  Kalman  filter  which  utilizes  the 
total  state  space  formulation.  The  filter  Is  mechanized  in  a  feedback 
configuration  (as  opposed  to  feedforward)  which  keeps  the  IKS  errors 
small  by  continually  feeding  back  a  correcting  signal  (see  Figure  1). 
Since  the  IKS  is  corrected  after  each  measurement  sample,  the  predicted 
error  states  and  measurement  differences  at  the  next  sample  tire  will 
be  zero.  Consequently,  there  is  no  need  to  propagate  the  error  state 
variable  estimates. 

A  discussion  of  the  relative  merits  of  both  direct  versus  Indirect 
and  feedforward  versus  feedback  configurations  Is  found  In  Reference  4 
pages  1-10.  It  is  not  the  Intent  of  this  report  to  dwell  on  the  develop¬ 
ment  of  Kalman  filter  theory,  rather  the  purpose  here  is  to  apply  this 
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Corrected 

Navigation 

Information 


Fig.  1.  Block  Diagram  of  Filter  in  Feedback  Configuration 
theory  to  a  practical  engineering  problem.  Therefore,  the  reader  who 
ia  either  unfamiliar  with  Kalman  filter  theory  and  its  applications 
to  inertial  navigation  systems  or  else  seeks  more  detailed  information 
on  this  subject  is  referred  to  References  2  through  5  in  the  bibliography. 

Limitations  of  the  Study 

It  is  realized  that  no  matter  vhat  system  reference  model  or  "truth 
model"  is  selected  it  will  only  be  an  approximation  since  the  complex 
real  world  dynamics  defy  exact  mathematical  description.  The  validity 
of  the  filter  design  is  sensitive  to  erroneous  reference  system  models 
and  statistics;  thus  great  care  must  be  taken  to  model  the  reference 
system  as  accurately  as  possible  in  the  computer  simulation  to  yield 
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XI •  Coordinate  Systems.  Transformation  Matrices,  and  Angular  Rates 

Coordinate  Frame  Definitions 

Three  different  coordinate  frames  are  defined  for  use  in  the  com- 
puter  simulation.  These  are  the  inertial  reference  frame,  the  Earth- 
fixed  frame,  and  the  navigation  frame.  The  following  paragraphs  will 
summarize  the  definitions  of  these  coordinate  systems. 

Inertial  Reference  Frame.  The  inertial  coordinate  system  is  fixed 
et  the  center  of  the  Earth  and  maintains  a  constant  orientation  with 
respect  to  inertlel  space.  The  coordinate  system  thus  defined  is  not 
trul>  an  inert lally-flxed  frame,  but  for  the  error  analysis  of  a  vehicle 
moving  near  the  Earth,  the  errors  introduced  by  this  definition  of 
inertial  space  are  negligible.  The  axes  of  the  coordinate  system  form 
an  orthogonal  right-handed  triad  with  the  x-axis  pointing  from  the 
center  of  the  Earth  through  the  Forth  Pole  in  alignment  with  the  Earth's 
spin  axis. 

Earth- fixed  Frame.  This  coordinate  system  is  identical  to  the 
inertial  reference  system  with  the  exception  thet  it  is  allowed  to 
rotete  with  the  Earth.  The  exes  of  this  coordinate  system  are  oriented 
vlth  the  x-axis  directed  outward  through  the  Korth  Pole,  the  y-axis 
directed  outverd  through  the  intersection  of  the  90  degree  West  meridian 
and  the  equator  (0,-90),  and  the  x-axis  directed  through  the  intersection 
of  the  Greenwich  meridian  and  the  equator  (0,0).  Since  this  frame  ro¬ 
tates  with  the  Earth,  any  point  on  the  Earth's  surface  can  be  specified 
in  terms  of  a  set  of  Euler  angle  rotations  in  the  Earth-fixed  frame,  A 
about  the  x-axis  (longitude)  and  0  about  the  y-exis  (latitude)* 

Navigation  Frame.  The  nevlgatloa  frame  is  the  coordinate  system  in 
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which  the  navigation  problem  la  aolved.  The  position,  velocity,  and 
attltuda  errors  are  expressed  in  navigation  coordinates.  Sines  the 
problem  involves  a  glab ailed  platform  (as  opposed  to  strapdovn)  which 
remains  essentially  locally  level  at  all  times,  the  z-axls  of  the  nav¬ 
igation  frame  remains  perpendicular  to  the  Earth  (positive  upward). 

There  are  several  choices  available  for  mechanization  of  the  locally 
level  system  in  terms  of  the  azimuth  angular  rata.  Among  these  are 
azimuth  wander,  constant  azimuth,  unipolar,  and  free  azimuth.  There 
are  certain  advantages  to  each,  for  example,  azimuth  wander  would  ba 
I’jed  in  navigation  near  the  polar  regions.  For  this  simulation,  however, 
constant  azimuth  mechanization  with  the  x-axis  of  the  system  always 
pointing  North  was  selected.  The  y-axis  is  than  chosen  as  pointing 
West  to  complete  the  right-handed  orthogonal  set. 

Transformation  Matrices 

In  order  to  transform  vector  quantities  such  as  positions  and 
velocities  from  one  of  the  above  coordinate  frames  to  another  it  is 
convenient  to  first  derive  a  sat  of  time  varying  transformation  matrices 
which  will  be  used  in  subsequent  development  of  the  simulation  aquations. 

Inertial  to  Earth  Transformation.  The  transformation  from  the 
inertial  reference  frame  to  the  Earth  fixed  frame  involves  simply  a 
rotation  about  the  x-axis  through  an  angle  equal  to  the  Earth  angular 
rate  (which  is  constant)  multiplied  by  the  time.  Projecting  vector 
components  of  the  Inertial  frame  ale  r  the  axes  of  the  Earth  frame  and 
installing  these  as  columns  of  the  t  3 formation  matrix  yields 
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Earth  to  Navigation  Trana  format  ion.  This  transfomation  matrix  is 
darlved  similarly  to  ths  previous  one  only  this  time  involving  two 
rotations.  The  first  one  through  an  angle  A  about  the  x-axls ,  and  the 
second  a  rotation  about  the  y-axis  through  an  angle  6. 


cose  0  -sinO 
0  10 
sine  0  cos6 


10  0 
0  cosA  slnA 
0  -sin A  cosA 


cos6  •sine  slnA  hcosA  sln6] 
1  cosA  1 


I 


slnA 


(2) 


[sine  |-slnA  cos6|  cose  cosA] 

Inertial  to  Navigation  Transformation.  This  transformation  ia 
simply  a  product  of  the  above  two. 

Ci®  -  C.“Cl®  O) 

Also  note  that  C.*  ■  (C.*)”1  »  (C.*)T  and  similarly  for  cA  C  *. 

V  X  X  HU 


Angular  Rates 

The  angular  rates  are  now  obtained  and  coordinatised  in  the  navi¬ 
gation  frame,  since  this  is  the  frame  in  which  the  computations  will  be 
performed.  A  profile  generating  computer  program  which  will  be  dis¬ 
cussed  in  detail  in  Chapter  VIII  is  used  to  calculate  and  store  values 
related  to  the  dynamics  of  the  aircraft  and  satellites  at  each  time 
increment  of  numerical  integration.  The  velocity  and  acceleration 
components  (as  well  as  other  pertinent  values)  are  stored  and  available 
in  navigation  coordinates  and  are  thus  used  in  the  following  development 
of  anguler  rates.  These  angular  rates  and  ecceleratlons  will  be  re¬ 
quired  in  the  propagation  of  the  Plnaon  inertial  system  error  model 
plant  states.  Also  note  that  the  aircraft's  distance  from  the  Earth's 
center  is  approximated  as  being  equal  to  the  radius  of  the  Earth.  This 
is  Justified  by  the  fact  that  the  altitude  of  the  aircraft  is  on  the 
order  of  two  miles  as  compared  to  approximately  4000  miles  for  the 
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radius  of  ths  Earth. 

The  angular  rate  of  the  Earth  vlth  respect  to  the  inertial  frane 
and  coordinatlzed  in  the  navigation  frame  is 
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To  get  the  angular  rate  of  the  navigation  frame  with  respect  to 
the  Earth-fixed  frame  using  the  stored  values  of  aircraft  velocity 
coordinatlzed  in  the  navigation  frame 

(1)  rotate  through  an  angle  A  about  Ex 

(2)  rotate  through  an  angle  6  about  Ey 
where  Ex  and  Ey  are  defined  in  Figure  2.  Then 

s«  *  *  *x  *  5  V  <5) 


and 


R  cose 


(6) 


and 


(7) 


where  ex  and  Cy«  are  unit  vectors  along  the  axes  about  which  rotation 
occurs.  The  angular  rate  of  the  navigation  frame  with  respect  to  the 
Earth  in  navigation  coordinates  becomes 


px  ■  A  cos  0  ■  -  Vy/R 
Py  -  0  -  Vx/R 

PE  ■  A  sin  0  ■  -  Vy/R  tan  0 


(8) 
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where  R  Is  the  radius  of  the  Earth.  The  angular  accelerations  are  now 
obtained  by  taking  the  time  derivatives  of  the  angular  rates. 

m 

P, 


•n 

“en 


px  “  (-  VR)  -  -  VR 


py  ■  it  (VR>  -  VR 

.  d  V-  iL,  V_  ,  , 

pj  ■  ^  (-  tan  0)  ■  -  tan  0  -  ■jp  0  sec2  6 


(12) 

(13) 

(14) 

(15) 


The  angular  rate  of  the  navigation  frame  with  respect  to  the  inertial 
frame  expressed  in  navigation  coordinates  is  then. 
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(Reference  6:3-27  thru  3-43) 
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III.  Satellite  Geometry  and  Range  Measurement  Equation 

In  this  chapter ,  theory  which  la  peculiar  to  the  use  of  Navigation 
satellites  as  updating  aids  to  the  INS  is  developed.  First,  the  range 
neaaurement  equation  is  derived.  Next,  a  method  for  determining  whether 
or  not  a  satellite  is  observable  or  “in  view”  is  presented.  Finally, 
the  satellite  motion  generator,  a  computer  routine  which  calculates 
satellite  orbital  elements  and  unit  vectors  from  the  aircraft  to  the 
satellite  as  functions  of  time,  is  discussed. 

Range  Measurement  Equation 

The  range  divergence  equations,  characterising  the  range  measuring 
process,  are  generated  by  the  user  from  a  combination  of  INS  and  satellite 
information.  The  range  measurement  process  involves  the  comparison  of 
a  measured  value  of  range  against  a  predicted  value  of  range.  The 
measured  range  to  a  satellite  is  determined  by  measuring  the  incremental 
phase  shift  between  the  user  and  satellite  clocks  which  were  synchronised 
at  an  earlier  time.  The  computed  range  is  obtained  from  satellite  ephe- 
meris  data  and  user  INS  supplied  position  information.  Both  the  measured 
and  computed  range  values  contain  errors;  however,  if  the  measured  and 
computed  values  are  subtracted,  the  difference  will  contain  only  the 
errors.  This  difference  is  called  the  range  divergence.  If  these  errors 
can  be  modeled  as  the  outputs  of  linear  systems  driven  by  white  Gaussian 
noise,  then  a  Kalman  filtar  can  be  constructed  to  estimate  these  errors 
and  greatly  Improve  upon  the  accuracy  of  the  raw  range  data. 

In  order  to  avoid  the  notational  Inconvenience  of  using  subscripts 
to  keep  track  of  which  satellite  is  being  referred  to,  consider  the  case 
where  a  single  satellite  is  observed.  The  results  are  identical  for  any 
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on*  of  the  27  satellites  and  therefore  easily  extended.  The  range  vector 
of  Interest  Is  the  vector  r  from  the  aircraft  (or  user)  to  the  satellite. 
It  Is  related  to  the  two  vectors  r_  and  r„  which  are  illustrated  and 

— ■  — O 

defined  in  the  following  figure. 


I 


jr  -  position  vector  from  user  to 
satellite 

rg  ■  position  vector  from  Earth's 
center  to  aircraft 

rg  ■  position  vector  from  Earth's 
center  to  satellite 


Figure  3.  Definition  of  Range  Vector 

It  follows  that 

£  -  la  -  la  U7) 

r's  111  -  lla  -  I«l  <18> 

r  “  \/i  ’  £  -  \/(ls  "  I a>  *  <Is  ~  £*)  <19> 

Tue  measured  range  to  the  satellite,  r',  is  composed  of  two  parts. 


r'  -  r  +  fir'  (20) 

where  r  is  tne  true  range,  and  fir*  is  the  error  in  the  measured  range. 
The  computed  range  to  the  satellite,  r*,  is  also  composed  of  two 


parts , 


r*  ■  r  +  fir* 


(21) 


where  r  is  the  true  range,  and  fir*  is  the  error  in  the  value  of  the 
computed  range.  The  quantity  that  is  being  observed  is  then  the 


13 


GA/EL/74K-1 

difference  of  these  two  range  values  and  is  called  the  range  divergence, 
Ar. 

Ar  ■  r'  -  r*  ■  6r'  -  fir*  (22) 

Errors  in  Computed  Range.  The  computed  satellite  position  is  in 
error  due  to  errors  in  the  ephemeris  data,  while  the  IN'S  errors  account 
for  tne  uncertainty  in  the  Aircraft's  position. 

I«  "  J*  +  fir*  (23) 

r*  ■  r^  +  fir£  (24) 

The  error  equation  is  then  obtained  by  noting  from  equation  (19) 

that 

(r*)2  -  r  *  •  r  *  (25) 

Taking  the  differential  of  both  aides  of  the  equation, 


2r*fir*  ■  r*  •  fir*  +  fir*  •  r* 


(26) 


or. 


fir*  -  *  fi£*)  '  (  ~*£*)  *  fi£* 


(27) 


Noting  that  the  quantity  in  parentheses  on  the  right  hand  side  of  the 
above  equation  is  simply  a  unit  vector  from  the  aircraft  to  the  satellite 


i  =  —  r* 
l“* 


(28) 


and  also  that 


fir*  ■  fir*  -  ;-r* 
—  — s  -a 


(29) 


the  computed  error  is  then 


fir*  ■  ir*  ♦  [ 


(30) 


Satellite  orbital  parameters  are  undated  by  the  ground  tracking  network 
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on  a  periodic  basis  and  relayed  to  the  user  along  with  the  range  data. 
Thla  ephemeris  data  is  highly  accurate;  any  uncertainties  in  computed 
satellite  range  can  be  accounted  for  by  increasing  the  satellite  clock 
phase  error.  Thus,  it  is  assumed  that  is  approximately  zero,  and 
the  computed  range  error  is 


6r*  *  -  ir*  •  dr*  (31) 

The  computation  of  the  above  equation  requires  values  for  the  unit  vector 

from  the  aircraft  to  the  satellite  and  current  values  for  the  x,  y,  and 

T 

z  INS  position  error  states  (6r*  ■  (Ax,  Ay,  At]  )•  Actually,  the  root- 
mean-squared  (RMS)  values  of  the  covariance  of  the  position  errors  are 
used,  since  this  is  a  stochastic  process  simulation;  this  statistical 
aspect  of  the  problem  will  be  discussed  in  later  chapters. 

Errors  in  Measured  Range.  Modeling  of  the  range  measurement  error 
requires  a  knowledge  of  the  various  error  sources  which  are  contained  in 
the  measurement  and  fitting  these  error  sources  with  empirical  data. 

The  model  used  in  this  simulation  is  a  somewhat  simplified  version  of 
the  one  found  in  Reference  1.  It  is  a  linear  combination  of  three 
components  for  eacn  satellite  measurement  corrupted  by  white  Gaussian 
noise  (v).  Each  of  the  separate  components  is  a  linear  system  driven 
by  white  Gaussian  noise.  These  linear  models  will  be  discussed  in  detail 
in  Chapter  V. 

The  range  measurement  error  is  modeled  as 

dr'  •  c6Tu  -  cdTg  +  db  +  v  (32) 

where 
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c  ■  the  speed  of  light 

6TU  ■  user  clock  phase  error 

6TS  ■  satellite  clock,  phase  error 

6b  ■  range  bias 

v  *  measurement  noise 

The  bias  term  In  the  above  equation  accounts  for  the  minor  effect 
of  both  tropospheric  delay  and  velocity  of  light  bias  uncertainties  in 
each  of  the  4  satellite  range  measurements.  The  error  due  to  ionosphere 
delay  is  a  function  of  the  elevation  angle  and  on  the  order  of  15  feet. 

A 8  will  be  explained  in  Chapter  V,  this  effect  is  assumed  to  be  in¬ 
cluded  in  the  satellite  clock  phase/range  error  (6TS). 

Substitution  of  equations  (31)  and  (32)  into  equation  (22)  yields 
the  final  form  of  the  range  divergence  equation. 

Ar  •  ir*  •  fir*  +  cfiT,  -  cfiT  +  fib  +  v  (33) 

A  minimum  of  four  range  divergence  equations,  or  in  other  words 
at  least  four  satellites,  are  required  as  observables  to  correct  for 
the  three  components  of  position  and  the  clock  phase  or  tine  difference. 
Also  note  that  for  notatlonal  convenience  the  asterisk  (*)  will  be 
dropped  in  subsequent  equation  development. 

Satellite  In- View  Criterion 

In  order  to  obtain  measurements  from  a  given  satellite,  the  sate¬ 
llite  must  be  observable  by  the  user;  that  is,  it  must  be  at  some 
specified  minimum  angle  above  the  aircraft's  horiton  for  a  useful  signal. 
This  mlnlnun  angle  is  arbitrary  and  dependent  upon  the  capabilities  of 
the  user  equipment.  A  nominal  value  of  ten  degrees  was  selected  for 

I 

implementation  in  this  study.  This  ln-vlew  criterion  along  with  the 
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selection  of  27  for  the  number  of  satellites  globally  deployed  insures 
that,  regardless  of  the  user's  location,  a  reasonable  number  (say  seven 
or  eight)  satellites  will  always  be  observable,  from  which  a  "best  set" 
of  the  required  four  satellites  may  be  chosen. 

A  method  for  determining  whether  or  not  a  satellite  is  in-view  using 
information  calculated  in  the  satellite  orbit  generating  computer  program 
is  now  presented. 


I 


-  Satellite 


Figure  4.  In-View  Criterion  Geometry 
■  minimum  angle  of  elevation  for  useful  satellite  signal 

■W  ■  «»*  -  (*> 

r,  ra,  ra,  and  Ca  are  as  defined  in  the  previous  section. 
ra  is  most  readily  expressed  ir  the  navigation  frame. 


where  R  ■  radius  of  Earth,  h  -  aircraft  altitude,  and  the  superscript 
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denotes  the  frame  in  which  the  vector  is  coordinetized.  Since  re  is 
derived  from  the  ground  track  latitude  and  longitude  of  the  satellite 
in  the  orbit  generator  and  readily  available,  the  vector  r  from  the 
aircraft  to  the  satellite  coordlnatized  in  the  navigation  frame  is 
written  as 

nnn  n  e  _ 

£  •  Is  I«  " 

A  unit  vector  along  £n  is  given  by 


(37) 

(33) 


Noting  that  the  z-component  of  the  above  unit  vector  is  simply  the  sine 
of  the  elevation  angle  or  the  cosine  of  D, 


It  follows  that  if  the  unit  vectors  from  the  user  to  the  satellite 


expressed  in  navigation  coordinates  are  computed,  the  in-view  criterion 
becomes 


cos  D  >  cos  Dmax  (41) 

If  ir^  >  cos  Dnvax  satellite  is  in  view. 

If  ir  <_  cos  Dmiy  satellite  is  not  in  view. 

For  the  arbitrarily  selected  minimum  elevation  angle  of  ten  degrees 
this  criterion  requires  the  z-component  of  this  unit  vector  to  be 
greater  than  cos  80*. 
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ir  >  cos  80*  (cos  80*  =  0.17A)  (A2) 

X 

A  sample  output  from  the  satellite  orbit  generator  program  listing 
satellites  which  are  in  view  at  a  particular  time  Instant  is  Included 
at  the  end  of  this  chapter. 

Satellite  lotion  Generator 

It  is  apparent  from  previous  sections  that  in  order  to  simulate  an 
INI  flight,  a  satellite  motion  generator  Is  required  to  provide  necessary 
information  regarding  the  satellites'  orbital  elements.  This  section 
presents  the  equations  necessary  to  calculate  required  parameters  of  the 
satellites  and  also  describes  the  selected  arbitrary  initial  conditions 
on  the  satellite  constellations.  It  should  be  noted  that  this  particular 
configuration  is  a  leading  candidate  for  implementation  in  the  proposed 
Hi  I  system. 

Total  deployment  consists  of  three  rings  or  "constellations"  of 
nine  satellites  each.  All  satellite  orbits  are  assumed  to  be  circular; 
in  fact,  the  orbital  period  (and  thus  the  orbital  velocity  and  altitude) 
of  all  satellites  is  assumed  constant  and  equivalent.  Since  global 
coverage  is  desired,  the  satellites  on  any  given  ring  are  equally  spaced; 
thus,  the  circular  arc  between  any  two  adjacent  satellites  on  a  ring 
subtends  a  central  angle  of  forty  degrees.  A  satellite  is  identified 
by  a  two  digit  code,  the  first  digit  (1  through  9)  indicates  which 
satellite  of  the  nine  on  the  ring  is  referred  to,  and  the  second  digit 
indicates  the  particular  constellation  (1,  2,  or  3).  Thus,  satellite 
62  is  the  sixth  satellite  on  the  second  ring. 

The  Euler  Angles.  To  specify  the  orientation  of  any  one  satellite 
with  respect  to  the  Earth-fixed  frame  requires  three  parameters;  the 
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most  convenient  parameters  are  the  Euler  angles,  from  which  the  direc¬ 
tion  cosines  or  unit  vector  to  the  satellite  may  be  determined.  Figure 
5  shows  two  unit  orthogonal  right-handed  triads  (i,  p  k)  and  (1^,  K) 
with  origins  at  point  0.  The  triad  (i,  k)  has  Jc  pointing  to  the 
satellite,  and  the  triad  (I,  J_,  K)  represents  the  Earth-fixed  frame. 

The  first  Fuler  angle  C  is  the  angle  between  jL  and  1^  and  is 
the  angle  of  inclination.  The  second  angle  n  is  the  angle  between  the 
plane  (i,  p  and  the  plane  U,  J) ;  this  is  the  angle  resulting  from 
the  Earth's  rotation.  The  third  angle  c  is  the  angle  between  the  plane 
(i,  3)  an^  the  plane  (I.,  p;  this  is  the  angle  resulting  from  the 
satellites  trotlon  about  the  Earth. 

Initially  let  the  triad  (1,,  k)  coincide  with  (I,  £,  K) .  The 
triad  (p  p  k)  can  be  brought  to  the  general  position  shown  in  Figure 
5  by  applying  the  following  rotations  in  order  i 

(i)  A  rotation  n  about  p  this  brings  the  movable  triad  (i,  p  V) 
into  coincidence  with  (I,  ' ,  K' ) . 

(ii)  A  rotation  c  about  K' ;  this  brings  (i,  p  je)  into  coincidence 
with  (i,  J",  U'K 

(ill)  A  rotation  c  about  1,;  this  brings  (_i,  p  k)  into  the  required 
final  position. 

It  is  observed  that  all  possible  positions  of  the  body  can  be  ob¬ 
tained  by  assigning  values  to  C,  n,  S  in  the  ranges 

0  <_  £  <,  w  0  ^  n  <  2n  0  <,  C  <  2w 

Let  A  be  the  rotation  matrix  trn  ;;  forming  a  unit  vector  from  the 
general  coordinate  frame  (jL,  p  irth-fixed  frame.  Then, 
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a  . 
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*  « 

•  m 

31 

32 

33 

m 

(43) 


Where  I,  J,  U,  i,  j,  k  represent  the  scalar  components  along  their 
respective  axes.  Since  the  unit  vector  to  the  satellite  is  coincident 
with  k  (i.e.,  i  ■  j  •  0,  k  •  1)  this  simplifies  to. 


(44) 


These  direction  cosine?  may  be  obtained  by  resolving  vectors  In  Figure 
5.  Alternately,  the  cosine  lav  of  spherical  triangles  may  be  applied 
at  this  point  (See  r.ef  7:17-1?).  Using  either  method  the  following 
values  are  obtained. 


al3  •  sin  C  sin  C  (45) 

a23  a  -(sin  n  cos  x.  +  cos  n  sin  C  cos  E)  (46) 

a33  ■  cos  n  cos  c  -  sin  n  sin  Q  cos  E  (47) 

(Ref  8:259-261). 

*13*  *23*  *n<*  *33  rePresent  thG  components  of  the  unit  vector  along 
r*«  the  vector  from  the  center  of  the  earth  to  the  satellite  expressed 
in  Earth  coordinates. 

Initialisation  of  Satellites.  Tis  initial  conditions  and  constant 
parameters  of  the  satellite  orbits  .  i?  given  In  the  following  two  tables. 
Note  tuat  m  refers  to  the  mth  satci  te  on  ti  e  designated  ring;  for 
example,  the  entry  of  ml  in  the  r  initial  conditions  represents 

the  remaining  eight  satellites  o*  c...  .xist  constellation.  The  numbers 
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which  are  missing  from  the  table  (represented  by  — )  are  dependent  upon 
the  Initial  values  of  5,  n»  and  (  which  are  specified.  The  latitude 
and  longitude  values  refer  to  the  ground  track  of  the  satellites. 

TABLE  1 

Orbital  Design  Constants 

Orbital  Period  8  Hrs.  (all  satellites) 

Angle  of  Inclination  (?)  55*  (all  three  satellite  rings) 

Altitude  7496  n.  mi.  (all  satellites) 

TABLE  II 

Navigation  Satellite  Initial  Conditions 


Initial 

Latitude 

(deg) 

Initial 

Longitude 

(deg) 

• 

?o 

(deg) 

11 

0 

0 

0 

ml 

— 

— 

-40(m-l) 

12 

0 

120 

0 

m2 

— 

— 

-40(m-l) 

13 

0 

-120 

0 

m3 

_ 

-40(m-l) 

These  satellites  may  be  initialised  at  any  given  time  relative  to 
the  flight  time  since  they  are  in  reality  not  yet  in  orbit  and  their 
orbits  remain  arbitrary.  For  computational  ease,  the  satellites  will 
be  given  the  above  initial  conditions  at  the  start  of  the  INI  flight 
(t  ■  0).  however,  if  the  satellites  were  actually  in  orbit,  then 
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initial  positions  would  be  dependent  upon  the  user  take-off  time.  A 
pictorial  representation  of  the  satellite  orbits  is  shown  in  Figure  6. 

Computational  Sequence.  The  sequence  of  equations  required  to  com¬ 
pute  unit  vectors  for  each  of  the  27  jrn  is  now  presented. 


K  -  55* 

n  -  120* (n-1)  -  (l*/240  sec) t 
C  -  -40*(ra-l)  +  Ct 


(48) 

(49) 

(50) 


Where  n  ■  satellite  designator,  n  ■  ring  designator,  and  l*/240  sec 
represents  the  rotational  speed  of  the  Earth.  Also, 


C  ■  /G#/rae  -  orbital  speed  '  the  satellite  (51) 

The  components  of  r?  are  now  computed  us^ng  the  Eulerian  Angles. 


(52) 


*13 

sin  C  sin  £ 

*23 

m 

-(sin  n  cos  (  +  cos  n  sin  c  cos  O 

*33 

cos  r\  cos  c  -  sin  n  ®in  £  cos  5 

■  m 

•  a 

The  ground  track  latitudes  and  longitudes  arc  given  by 


8  -  tan"1  (a13/Aa23)z  +  (a^)^  ) 


A  -  tan"1  (-a23/a33) 


(53) 

(54) 


The  required  components  of  the  unit  vector  along  rn  nay  now  be  computed. 

(55) 

(56) 


rj 


Cn  re 

e 


rn  - 


r  -  /(rx)z  +  (r  )*  +  (r£)2 


(57) 
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Navigation  Satellite  Orbits 


ir  >  cos  80*  ?  (59) 

*■* 

The  three  components  of  the  unit  vector,  the  latitude,  longitude, 
and  in-view  criteria  for  each  of  the  27  satellites  is  calculated  and 
stored  at  each  tine  increment  of  numerical  integration  along  the  flight 
path.  A  sample  output  of  these  values  is  presented  in  Table  III.  This 
table  shows  the  output  of  the  profile  generator  program  for  the  27 
satellites.  The  first  block  is  ring  1  consisting  of  nine  satellites. 
Columns  1,  2,  and  3  are  the  x,  y,  and  z  components  of  the  unit  vector 
from  the  aircraft  to  the  satellite  expressed  in  the  navigation  frame. 
Columns  4  and  5  are  the  satellite  ground  track  latitude  and  longitude 
respectively  (in  degrees).  Column  6  is  the  in-view  criterion,  the 
component  of  the  unit  vector  in  the  z  direction.  If  this  value  is 
greater  than  cos  80*  5  0.17,  then  column  7  will  contain  a  T  (for  True) 
Indicating  the  satellite  is  in-view.  The  2nd  and  3rd  blocks  represent 
rings  2  and  3.  Note  that  the  satellites  in  view  of  this  particular 
time  instant  are  satellites  52,  62,  72,  13,  83  and  93. 
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Table  III 


i^5i\TcopV: 


Satellite  Position  and  In-View  Criterion  at  Tire  -  360  sec 


Satellite 

Number 


Ring  1 


11 

21 

31 

41 

51 

61 

71 

81 

91 


.376 

-.996 

u44 

3.7 

1.1 

-.34 

F 

.714 

-  •  b95 

-  *  07  8 

-26.4 

-22.7 

-.  18 

F 

.9  73 

u  77 

-.217 

-52.5 

-87.2 

-.22 

F 

•  766 

.517 

-.382 

-47.7 

-121.2 

-.3  8 

F 

•  243 

•  833 

-.495 

-15.9 

-186  .e 

-.5  0 

F 

•  363 

•  773 

-.521 

12.6 

165.5 

-.52 

F 

.817 

•  366 

-.447 

42.5 

i?e.7 

-.45 

F 

.923 

-.243 

-.299 

54.6 

75.0 

-.3  3 

F 

•  5  87 

-.7  93 

-.137 

3b.  1 

27.9 

-.14 

F 

Ping  2 


12 

22 

32 

42 

52 

62 

72 

82 

92 


13 

23 

33 

43 

53 

63 

73 

83 

93 


•  447 

•  U5 

•  360 
■•724 
•.836 
•*  466 

.203 

.740 

.723 


■*368 

•  5  44 
••448 
••  2C6 

•  1C? 

•  419 

•  625 

•  5  EC 

•  101 


196 

-.873 

3.7 

121.1 

-.87 

F 

129 

-.991 

-28.4 

56.3 

-.95 

F 

4o7 

—  .633 

-52.5 

52.8 

-.83 

F 

561 

-.403 

-47.7 

-11.2 

-.40 

F 

486 

.2  56 

-19.5 

-46.8 

•  26 

T 

458 

.  6  83 

12.6 

-70.5 

•  ee 

T 

496 

•  6  14 

42.5 

-101.3 

.81 

T 

656 

•  145 

54.6 

-161.0 

•  16 

F 

493 

-.483 

35.1 

147.9 

-.48 

F 

3 

645 

.  36  Li 

3.7 

-116.9 

•  36 

T 

780 

-.  211 

-28.4 

-143.7 

-.31 

F 

442 

-.777 

-52.5 

172.8 

-.78 

F 

007 

-.578 

-47.7 

108.8 

-.98 

F 

413 

-.504 

-19.9 

73.2 

-.9? 

F 

716 

-.558 

12.6 

49*5 

-.56 

F 

772 

.  033 

42.5 

18.7 

.03 

F 

351 

•  7 1  >■ 

54.6 

-41.0 

•  71 

T 

391 

•  9  lb 

35.1 

-92.1 

•  91 

T 

O 
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IV.  Kalman  Filter  Equations 

Design  of  an  INI  Kalman  filter  requires  extensive  computer  simula¬ 
tion.  This  chapter  Is  a  discussion  of  the  equations  which  are  required 
not  only  for  the  mechanization  of  the  filter  but  also  those  which  are 
necessary  to  simulate  the  dynamics  of  the  user,  or  aircraft,  and  the 
driving  error  sources. 

The  method  of  covariance  analysis  will  be  the  principal  tool  used 
In  solution  of  the  problem.  In  this  analysis,  arbitrary  Initial  condi¬ 
tions  on  the  diagonal  elements  of  the  covariance  matrix,  F,  are  specified 
and  the  off-diagonal  terms  are  assumed  to  be  zero  Initially.  The 
covariance  is  a  measure  of  the  uncertainty  in  the  knowledge  of  the  true 
values  of  the  components  of  the  state  vector.  In  the  simulation,  the 
covariance  matrix  of  both  the  system  and  the  filter  are  propagated  for¬ 
ward  in  tine  by  numerical  Integration  techniques.  When  the  specified 
update  time  is  reached,  the  best  estimates  of  the  states  are  determined, 
and  control  is  applied  to  the  system  to  adjust  the  values  of  the  state 
variable  to  the  best  estimate  obtained  with  the  Kalman  filter.  The 
square  root  of  the  individual  diagonal  elements  of  the  system  covariance 
matrix  are  then  plotted  as  a  function  of  time  to  yield  a  basis  for 
comparison  of  filters. 

In  this  simulation,  the  error  statistics  are  propagated;  i.e.,  the 
standard  deviation  of  the  noise  is  always  supplied  when  a  noise  value  is 
required.  This  is  in  contrast  to  a  "Monte  Carlo"  type  simulation  where 
an  actual  sequence  of  noise  values  is  measured  and  recorded  for  use  in 
the  simulation.  This  is  possible  because  the  covariance  is  Independent 

of  actual  measurement  values,  and  can  be  computed  without  generating  a 
sample  sequence  of  measurement  values. 
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System  Model  Equations 

The  basic  equations  used  in  this  process  are  the  differential  equa¬ 
tions  that  describe  how  the  inertial  navigator  errors  propagate  with 
time.  The  equations  are  formulated  into  a  set  of  first  order,  linear 
differential  equations,  driven  by  white  Gaussian  noise.  Linear  measure¬ 
ments  are  made  upon  the  actual  system  variables,  and  these  are  corrupted 
by  white  Gaussian  noise.  It  is  assumed  that  the  equations  representing 
a  detailed  model  of  the  By»tem  are  as  follows* 


where 


xff  ■  FgXg  +  G.u.  (60) 

Xg  is  an  nj  vector  denoting  the  true  state 
Fs  is  an  n1  x  n.  system  dynamics  matrix 
Gg  is  an  n^  x  m^  gain  matrix 

Ug  is  an  m}  vector  of  white  noise  Inputs  with  zero  mean  and 
variance 


E[u(i)u(j)T)  -  j  (61) 

where  1  and  j  are  instants  in  time. 

The  observations  obtained  from  external  references  can  be  described  by 
the  linear  measurement  vector  equation 


*.  -  H  x  +  v  (62) 

where 

Sg  is  an  r  vector  of  measurements 
Hg  ia  an  r  x  n{  measurement  matrix 

Vg  is  an  r  vector  of  white  noise  Inputs  with  zero  mean  and 
variance 
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E[v(i)v(j)T]  -  f*8™  1  "  J  (63) 

L  0  if*  J 

It  Is  assumed  that  the  system  noise  and  measurement  noise  are  uncorrelated 
for  all  time. 

E[u(i)v(j)T]  -  0  all  i,J  (64) 


Filter  Equations 

Tne  above  set  of  equations  are  assumed  to  be  a  complete  and  accurate 
mathematical  description  of  the  INI  system  dynamics  and  measurement  equa¬ 
tions  for  the  purpose  of  simulation.  It  Is  also  the  set  of  equations 
which  would  be  utilized  In  the  design  of  a  fully  optimal  Kalman  filter. 
However,  due  to  the  computational  burden  of  the  optimal  filter,  a  sub- 
optimal  filter  design  Is  obtained  by  reducing  the  dimension  of  the  state 
vector.  The  states  that  are  eliminated  are  the  ones  that  least  affect 
the  accuracy  of  the  mathematical  description  of  the  IIII.  This  suboptlnal 
filter  can  then  be  implemented  with  an  aircraft  on-board  computer. 

The  suboptimal  filter  structure  is  represented  as  follows  where 

A 

X£  Is  the  design  filter  state  and  xf  is  the  filter  best  estimate  of  the 
design  filter  state. 

if  "  Ff*f  +  Gf«f  <65> 

where 

Xf  is  an  n2  vector 

Ff  Is  an  n2  x  n filter  dym  .:iics  matrix 
Gf  is  an  n2  x  m2  gain  matri .. 

Uf  is  an  m  vector  of  v'nJt'  Inputs  with  zero  mean  and 

variance 
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*f  m*i  +  Kfl*.  -  H^]  (73) 

where 

A 

£jr  le  an  n2  vector  denoting  the  best  estimate 
Pf  is  the  filter  covariance  matrix 

-  superscript  indicates  the  time  instant  before  update 
+  superscript  indicates  the  time  instant  after  update 
Kf  Kalman  gain  matrix 
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The  filter  takes  the  actual  measurement  z  and  subtracts  from  it  the 

— s 

best  prediction  of  its  value  before  the  actual  measurement  is  taken, 

* 

H^Xf,  This  difference  is  then  passed  through  an  optinal  weighting  ma- 
trix  Kf,  and  used  to  correct  jcp  the  best  prediction  of  the  state  at  the 
time  instant  before  the  measurement  is  taken.  This  gives  the  best 
estimate  after  update.  This  estimate  in  propagated  to  the  time  of  the 
next  measurement  sample  with  the  above  equations  jk  and  P. 

These  recursive  relationships  are  initiated  from  the  assumed  Gaussian 
density  that  describes  the  apriori  knowledge  of  the  state. 


jc(0)  - 

(74) 

P(0)  -  P0 

(75) 

The  Kalman  filter  will  propagate  the  conditional  probability  density 
of  the  desired  quantities,  conditional  on  the  actual  measurements  taken. 

The  probability  density  function  of  a  Gaussian  noise  amplitude  takes  on 
the  shape  of  a  normal^  bell-shaped  curve.  This  assumption  of  Gaussian 
noise  amplitude  is  justified  by  the  fact  that  a  system  or  measurement 
noise  is  typically  caused  by  a  number  of  small  sources.  It  can  be 
shown  mathematically  that  when  a  number  of  random  variables  are  added 
together  the  sunned  effect  is  very  nearly  a  Gaussian  probability  density, 
regardless  of  the  shape  of  the  individual  densities.  Additionally,  the 
use  of  Causslan  densities  makes  the  mathematics  tractable.  The  first 
and  second  order  statistics  (mean  and  variance  or  standard  deviation) 
completely  determine  a  Gaussian  density.  Thus,  the  Kalman  filter,  which 
propagates  the  first  and  second  order  statistics,  Includes  all  informa¬ 
tion  contained  in  the  conditional  probability  density  (Ref  10:4-6  &  2:6-21). 

The  mean  or  expectation,  (p),  of  a  density  is  defined  as 
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09 

F^xj  *  p  -  /  xf(x)dx 


This  is  interpreted  as  the  weighted  average  of  the  values  of  X,  us  inf* 
the  probability  density  function  f(x)  a3  the  weighting  function.  This 
value  is  assured  to  be  zero  for  all  Gaussian  white  driving  noises  used 
in  this  simulation. 

The  variance  (a2)  of  a  density  is  defined  as 


Var^xj  -  a2  ■  /  (x-p)2f(x)dx 


o2  is  the  weighted  average  of  the  values  of  (x-p)2;  thus,  o2  is  a  mea¬ 
sure  of  the  spread  of  the  density.  (Ref  11:136-146).  It  is  a  direct 
measure  of  the  uncertainty:  the  larger  o  is,  the  broader  the  probability 
peak  is,  spreading  the  probability  weight  over  a  larger  range  of  x 
values.  For  a  Gaussian  density,  68,32  of  the  probability  weight  is 
contained  within  the  band  o  units  to  each  side  of  the  mean  (p)  which 
represents  the  area  under  the  normal  bell-shaped  curve  between  -o  and 
+o  and  95.4%  of  the  probability  weight  is  contained  between  -2a  and 


Equations  76  and  77  give  the  first  and  second  order  statistics  for 
the  scalar  case.  These  equations  arc  easily  extended  to  the  vector  case 
(as  required  in  this  study)  and  are 


e[x1  "  £  ■  /  *  *  *  ■  xf(x)dxj  *  *  *  dXjj 

L  J  —00  ju 


(76a) 


Cov[x]  -  P  -  /  *  *  •  /  (ii.  *  p)(x  -  jj)Tf (x)dxj  *  *  *  d xn  (77a) 
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V.  System  Model 

To  apply  the  Kalman  Filter  equations  developed  In  the  previous 
chapter,  a  reference  system  model  that  Is  a  good  approximation  to  the 
real  world  dynamics  Is  needed.  This  chapter  outlines  the  reference 
system  equations  selected  for  this  design  study.  The  chapter  Is  pre¬ 
sented  In  four  sections.  The  first  section  defines  the  44  error  states 
Incorporated  in  the  system  model  along  with  their  assumed  initial  condi¬ 
tions.  Also,  the  system  dynamic  matrix  Ffi  is  presented  In  partitioned 
matrix  form.  The  second  section  discusses  the  modeling  of  the  IMS  plant 
error  states.  The  third  section  presents  IMS  and  satellite  error  source 
models,  linear  systems  driven  by  white  Gaussian  noise.  The  last  section 
displays  the  measurement  equation  (62)  explicitly  In  matrix  form. 

State  Variable  Definition  &  Initial  Conditions 

Table  IV  presents  a  detailed  listing  of  the  44  states  utilized  in 
the  reference  system*  model.  The  Initial  conditions  on  the  IMS  error 
states  are  highly  arbitrary,  and  values  similar  to  those  used  in  other 
studies  are  selected.  The  Initial  conditions  on  the  accelerometer  and 
gyro  error  states  are  typical  of  an  inertial  navigation  system  in  the 
one  to  two  nautical  mile  per  hour  class.  The  initial  conditions  on  the 
user  clock  are  similar  to  those  in  Reference  1.  The  initial  conditions 
on  the  satellite  clock  model  are  classified;  therefore,  reasonable  order 
of  magnitude  numbers  were  selected  in  order  to  maintain  the  unclassified 
status  of  this  report.  It  should  be  pointed  out  that  whenever  an 
integrator  is  required  in  the  simulation  diagram,  the  output  of  this 
integrator  is  a  system  state  variable.  Thus,  the  addition  of  an  inte- 

I 

grator  to  an  error  model  increases  the  dimension  of  the  error  state 
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Table  IV 

System  44  State  Vector  Definition 


Error 

State 

Symbol 

Definition 

RMS  Initial 
Condition 

INS  Plant  Error  States 

1 

Ax 

x  position  error 

3000  ft 

2 

Ay 

y  position  error 

3000  ft 

3 

As 

s  position  (altitude)  error 

300  ft 

4 

Ax 

x  velocity  error 

2  ft/sec 

5 

Ay 

y  velocity  error 

2  ft/sec 

6 

As 

z  velocity  error 

0.5  ft/sec 

7 

♦x 

x  attitude  error 

0.14  mlllirad 

8 

s 

y  attitude  error 

0.14  mlllirad 

9 

♦. 

z  attitude  (heading)  error 

2.0  mlllirad 

INS  Error  Sources 

10 

°rx 

x  random  accelerometer 
noise — expotentially 
correlated,  t  ■  600  sec 

20  vg’s 

11 

°ry 

y  random  accelerometer 
noise  expotentially 
correlated,  t  ■  600  sec 

20  yg's 

12 

ars 

z  random  accelerometer 
noise — expotentially 
correlated,  t  ■  600  sec 

20  vg's 

13 

“bx 

x  acceleoaetcr  bias 

66  tig's 

14 

“by 

y  acceleroact  r  bias 

66  vg's 

15 

“bs 

z  acceleronec .. r  bias 

132  vg'a 

16 

*rx 

x  random  r* " '  '  t 

expotenti:’’  rlated, 

t  ■  3600  sec 

0.012  deg /hr 
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Table  IV 

System  44  State  Vector  Definition 


Error 

State 

Symbol 

Definition 

RMS  Initial 
Condition 

17 

cry 

y  random  gyro  drift 
expotentlally  correlated ( 
t  •  3600  sec 

0.012  deg/hr 

18 

cr* 

z  random  gyro  drift 
expotentlally  correlated, 

T  *  3600  sec 

0.012  deg/hr 

19 

ebx 

x  gyro  bias 

0.025  deg/hr 

20 

eby 

y  gyro  bias 

0.025  deg/hr 

21 

cbz 

s  gyro  bias 

0.04  deg/hr 

User  Clock  Errors 

22 

«TU 

uner  clock  phase/rsnge  error 

1C, 000  ft 

23 

*23 

user  clock  frequency  offset 

5  ft/sec 

24 

x2k 

long  term  stability  error  ' 

5  x  10“6  ft/sec2 

Satellite  Clock  Errors 

25 

6T 

4  si 

sat.  clock  #1  phase/range 
error 

100  ft 

26 

*26 

sat.  clock  #1  frequency 
offset 

0.01  ft/sec 

27 

*27 

sat.  clock  #1  stability 
term 

1  x  10~7  ft/sec2 

28 

*28 

sat.  clock  #1  expotentlally 
correlated  noise,  t  ■  10000 

1  x  10"3  ft/sec 
sec 

29 

6T 

*2 

sat.  clock  #2  phase/range 
error 

100  ft 

30 

X30 

sat.  clock  #2  frequency 
offset 

0.01  ft/sec 
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Table  IV 

Syatem  44  State  Vector  Definition 
(Cont.) 


Error 

State 

Synbol 

Definition 

RMS  Initial 
Condition 

31 

*31 

aat.  clock  $2  stability  tens 

1  x  10“7  ft/sec2 

32 

*32 

sat.  clock  # 2  expotentially  1  x  10”  3  ft/sec 

correlated  noise,  t  ■  10000  sec 

33 

61.3 

sat.  clock  #3  phase/range 
error 

100  ft 

34 

*3H 

sat.  clock  #3  frequency 
offset 

0.01  ft/sec 

35 

*35 

sat.  clock  #3  stability  tern 

1  x  10"7  ft/sec2 

36 

*36 

sat.  clock  #3  expotentially  1  x  10“  3  ft/sec 

correlated  noise,  t  ■  10000  sec 

37 

6I.H 

sat.  clock  #4  phase/range 
error 

100  ft 

38 

*38 

sat.  clock  #4  frequency  offset  0.01  ft/sec 

39 

*39 

sat.  clock  #4  stability  tern 

1  x  10”7  ft/sec2 

40 

*H0 

sat.  clock  J4  expotentially 
correlated  noise,  t  ■  10000 

1  x  10”3  ft/sec 

Satellite  Range  3ios  Errors 

41 

«bl 

sat.  #1  range  bias 

10  ft 

42 

6b2 

sat.  #2  range  bias 

10  ft 

43 

6b3 

sat.  #3  range  bias 

10  ft 

44 

«bH 

sat.  #4  range  bias 

10  ft 
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variable  vector  by  one.  Also,  x,  y,  and  z  are  navigation  frame  axes. 

The  initial  covariance  matrix  P(0)  has  now  been  completely  speci¬ 
fied.  Its  diagonal  elements  are  the  squared  values  of  the  RMS  initial 
conditions  given  in  Table  IV.  The  remaining  off-diagonal  elements  are 
assumed  to  be  zero  initially.  Propagation  of  the  linear  variance  equa¬ 
tion  (70)  requires  additional  knowledge  of  the  two  matrices  F  and  Q1 
where 

q*  .  GQGT  (78) 

The  F#  matrix  is  partitioned  as  follows 
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Plant  Error  States 

The  plant  error  state  equations  are  the  differential  equations 
describing  the  natural  unforced  dynamic  response  of  the  errors  in  the 
inertial  navigation  system.  This  includes  nine  states:  x,  y,  and  z 
position,  velocity,  and  attitude  in  the  navigation  frame.  Additional 
state  variables  are  usually  added  to  these  for  the  purpose  of  damping 
the  inherently  unstable  vertical  channel.  Also,  an  altimeter  measure¬ 
ment  may  be  added  to  the  measurement  equations.  However,  since  the 
navigation  satellites  provide  position  information  along  all  three  of 
the  x,  y,  and  z  axes,  inclusion  of  these  additional  states  and  measure¬ 
ments  to  control  the  altitude  divergence  is  not  absolutely  necessary. 

In  actual  practice,  an  altimeter  will  be  provided  on  the  aircraft  and 
may  be  used  as  a  back-up  for  NAVSAT  equipment  failure  or  for  temporary 
measurements  during  the  period  a  satellite  goes  out  of  view  and  a  new 
one  is  being  acquired. 

There  are  various  models  of  these  nine  IMS  plant  states  available 
for  implementation.  The  Pinson  error  model  was  selected  as  it  was  the 
model  used  in  the  SAlIl’S  and  profile  generating  computer  programs  which 
will  be  explained  in  a  later  chapter.  A  derivation  of  the  Pinson  error 
model  is  given  in  Chapter  4  of  Reference  9.  This  model,  which  makes 
up  the  Fjj  matrix  is  shown  in  Figure  7. 

Error  Source  Models 

The  error  propagation  equations  given  in  Chapter  IV  were  developed 
under  the  assumption  that  the  system  disturbances  (u(t))  are  not 
correlated  in  time.  The  estimation  of  disturbances  which  have  signifi¬ 
cant  time  correlation  is  done  by  means  of  "state  vector  augmentation." 
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That  is,  the  dimension  of  the  system  state  vector  in  increased  by  in¬ 
cluding  the  correlated  disturbances  as  well  as  descriptions  of  their 
dynamic  behavior  in  the  appropriate  rows  of  an  enlarged  F  matrix.  Be¬ 
cause  these  quantities  are  random,  their  behavior  cannot  be  described 
deterministically.  Instead,  they  are  modeled  as  state  variables  of  a 
fictitious  linear  dynamic  system  which  is  excited  or  driven  by  white 
noise.  This  model  serves  two  purposes;  it  provides  proper  autocorrelation 
characteristics  through  specification  of  the  linear  system  and  strength 
of  the  driving  noise,  and  the  random  nature  of  the  signal  follows  from 
the  random  excitation. 

The  correlated  system  disturbances  utilized  in  this  study  are  each 
modeled  by  a  combination  of  one  or  more  of  the  several  types  of  basic 
error  models  described  in  Figure  8. 

Specification  of  the  block  diagram  models  of  these  error  sources 
implies  the  structure  of  the  F  and  Q*  matrices. 

The  random  constant  or  bias  is  a  non-dynamic  quantity  meant  to 
model  a  constant  of  unknown  amplitude.  It  is  simulated  as  the  output 
of  an  integrator  which  has  no  input  but  has  a  random  initial  condition. 

Its  constant  nature  is  indicated  by  cho  fact  that  the  corresponding 
rows  of  the  F  and  Q'  matrices  contain  only  zeros. 

Random  walk,  which  derives  its  name  from  an  illustration  involving 
a  man  who  takes  fixed-length  steps  in  arbitrary  directions,  is  simulated 
by  passing  white  noise  through  an  iut  .* orator.  In  this  case,  the  row  of 
the  augmented  F  matrix  contains  only  ;  ros.  However,  the  corresponding 
row  and  column  of  Q*  is  nonzero,  the  \lcnent  at  their  intersection  is  q. 

Random  errors  which  exhibit  -  .  ■  f  to  time-growing  behavior  may 
be  described  by  the  random  ramp  or  the  random  parabola  (which  is  formed 
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Fig.  8.  Error  Source  Models  for  Random  Variables 
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by  adding  a  long-term  stability  state  to  the  random  ramp) .  These  models 
will  Introduce  non-zero  elements  into  the  F  matrix,  and  when  driven  by 
white  noise,  into  the  Q'  matrix. 

The  expotentially  correlated  random  variable  whose  autocorrelation 
function  is  a  decreasing  expotential,  provides  a  reasonable  approximation 
for  a  band-limited  signal  whose  spectral  density  is  approximately  flat 
for  a  finite  bandwidth.  This  type  of  random  variable  introduces  a  single 
non-zero  diagonal  element  -8  into  the  F  matrix. 

8  -  ~  (80) 

where  t  -  correlation  time  constant  of  the  expotentially  decaying  func¬ 
tion. 

If  the  value  of  the  corresponding  element  in  the  initial  covariance 
matrix  P(0)  is  specified,  then  the  noise  value  required  to  drive  the 
fictitious  linear  system  for  computer  simulation  is  determined  as  follows. 
Since  q  is  constant,  in  steady  state  p  >  0  so  that  from  the  linear  vari¬ 
ance  equation 


0  -  FP  +  PFT  +  Q' 

(81) 

or,  in  the  scalar  case 

0  m  -gp  -  pfj  +  q 

(82) 

q  ■  2Bp 

(83) 

Therefore,  the  expotentially  correlated  random  variable  requires  a 
diagonal  element  in  Q'  whose  value  is  two  times  the  initial  covariance 
value  divided  by  the  correlation  time  (Ref  5:3-36  thru  3-49). 

INS  Error  Source  Models.  Through  extensive  testing  and  detailed 
knowledge  of  sensor  dynamics  many  imperfections  and  errors  of  inertial 
navigation  systems  are  removed  by  careful  design.  But  when  all  the  tests 
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for  predictable  errors  and  the  Ingenious  design  tricks  have  been  ex¬ 
hausted,  there  still  remain  errors  whose  source  defies  compensation.  The 
statistical  behavior  of  these  sensor  errors  can  be  obtained  through 
testing  and  fitting  curves  to  laboratory  experimental  data.  The  error 
models  chosen  for  the  gyros  and  the  accelerometers  are  typical  of  an  INS 
in  the  one  to  two  nautical  miles  per  hour  class.  They  are  both  modeled 
as  linear  combinations  of  a  random  bias  and  an  expotentlally  correlated 
random  variable  excited  by  white  noise.  Note  that  this  requires  two 
states  for  the  x  direction  gyro  drift  and  two  states  for  the  x  direction 
accelerometer  error,  and  similarly  for  the  y  and  z  directions.  There¬ 
fore,  modeling  of  the  INS  driving  errors  requires  12  additional  states 
In  the  system  error  state  vector.  Block  diagrams  and  corresponding 
elements  of  the  F  and  Q'  matrices  are  as  follows. 


I.C. 


u(t) 


a 


Fig.  9.  Accelerr::<  •  rr or  Source  Model 
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The  P  and  Q'  matrix  terms  are 


10  11  12  13-15 


Fig.  10.  Gyro  Error  Source  Model 
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The  F  and  Q*  matrix  terms  are  as  follows 


(90) 


<»16“2Bgpi6  (92) 

q16  -  2(1/3600  see)] (.012  deg/hr) (4.848  x  10“6  )]2 (93) 

q16  -  1.880  x  10"18  rad2/aec3  (94) 

Similarly , 

q17  "  ql8  "  q16  *95* 

Note  that  the  accelerometer  error  sources  are  additive  to  the  Ax,  Ay, 

and  Az  plant  state  differential  equations  requiring  the  non-zero  terms 
of  the  F12  matrix.  Similarly,  the  gyro  error  sources  ere  additive  to 
the  differential  equations  of  the  plant  attitude  error  states  $x»  4y» 
and  <|ix  requiring  the  non-zero  terms  of  the  FJ3  matrix. 

NAVSAT  Error  Sources.  Modeling  of  the  navigation  satellite  error 
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sources  requires  extensive  testing,  compiling  of  empirical  data,  and 
curve  fitting.  The  error  source  models  selected  for  this  design  study 
are  a  modified  version  of  those  found  in  Reference  1  and  include:  one 
user  clock,  four  satellite  clocks,  and  four  range  biases  (one  per  sate¬ 
llite). 

In  operation,  the  range  measuring  process  would  be  initialized  by 
synchronizing  the  user's  clock  (i.e.,  oscillator)  with  the  clock  signal 
received  from  the  satellite  at  some  arbitrary  starting  time.  Range 
Increments  would  then  be  measured  by  counting  incremental  phase  shift 
between  satellite  and  user  clocks  as  the  vehicle  moves.  If  either  clock 
drifts,  however,  erroneous  Incremental  phase  shifts  will  be  measured. 
During  the  tracking  process  the  satellite  clocks  are  in  effect  synchro¬ 
nised  to  the  master  tracking  station  clock  by  estimating  the  satellite 
clock  drift  model  coefficients.  These  coefficients  are  to  be  updated 
every  hour  and  relayed  via  the  satellite  user  link  to  the  user  for  use 
in  correcting  raw  incremental  range  measurements.  Thus  the  satellite 
clock  error  is  a  function  of  the  accuracy  of  the  drift  model  coefficient 
estimates  (Ref  1:2-14  thru  2-18). 

Reference  1  suggests  the  following  model  for  both  the  user  clock 
and  the  four  satellite  docks.  Note  that  although  the  structure  of  the 
user  and  satellite  clock  models  are  the  same,  the  user  clock  initial 
errors  are  orders  of  magnitude  greater  than  those  of  the  satellite  clock 
because  of  the  greater  accuracy  of  t.  ?  periodically  updated  satellite 
docks. 

The  mathematical  equation  repre.;  • .  ting  the  clock  error  is : 

6T  ■  CQ  +  ‘Jji.  2'.  -  +  e(t)  (96) 


47 


uA/ EE/7411-1 


where 

c(t)  +  fie  -  u(t)  (97) 

E[u(t) ]  -  0  (98) 

E[u2(t))  -  2eEUe(t)]2)  (99) 

The  block  diagram  of  the  fictitious  linear  system  simulating  this 
mathematical  equation  is  given  in  Figure  11. 


Fig.  11.  Clock  Error  Model 

The  values  suggested  for  the  satellite  clock  error  coefficients 
(C  ,  Ci,  C2,  c)  are  classified;  arbitrary  numbers  selected  for  the  de- 
sign  are  given  as  initial  conditions  in  Table  IV. 

The  contributions  to  the  system  F  and  Q'  matrices  by  the  satellite 
clocks  are  then 
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25  -  28 


29  -  32 


33  -  36 


37  -  40 


25-2G 


29-32 


55 


33-36 


37-40 


0  1  0  1  | 

0  0  1  0 

0  0  0  0, 

000  -0gl 

1 

0 

1 

0 

.0101 
o  |o  0  1  0 

0  0  0  0 
|°  0  0-6, 

0 

0 

0  10  1 

—  - 

0  |  o 

1 

0  0  10 

0  0  0  0 

0  0  0  -B. 

0 

0  I  0 

1 

0 

0  10  1 

0  0  10 

0  0  0  0 

0  0  0  -B 

q25"  20SP25  (101) 

^25  "  2(1/10000  «ec)(l  x  10-3  ft/*ac)2  (102) 

q25  "  2  x  10~10  *t2/««c3  (103) 

Similarly, 

q29  •  q33  •  q37  •  ^25  (104) 

Tha  suggested  user  clock  error  coefficient*  are  given  in  Table  V. 


Table  V 

Suggested  User  Clock  Error  Coefficients 


Error 

Coefficients 

RMS 

Value 

Correlation 

Time 

Co 

10000  ft 

- 

C| 

5  ft/sec 

- 

C2 

5  x  10“6  ft/sec2 

- 

• 

c 

2.5  ft/sac 

.02  sec 

(100) 


49 


GA/EE/74M-1 


These  numbers  are  typical  of  a  good  quality,  temperature  controlled 
quartz  crystal  oscillator  suitable  for  airborn  application.  These  values 
are  much  higher  than  the  satellites'  but  will  be  greatly  reduced  by  the 
inherent  calibration  process  of  the  i:>!X  Kalman  filter.  Random  drift  of 
the  user  clock  is  caused  predominantly  by  the  acceleration  sensitivity  of 
the  crystal  which  responds  to  aircraft  vibrations.  The  short  correla¬ 
tion  time  of  this  process  is  due  to  the  relatively  wideband  nature  of 
the  aircraft  vibration  spectrum  (Ref  1:2-17  thru  2-18). 

However,  the  correlation  time  of  the  c  state  variable  in  the  user 
clock  is  so  short  in  comparison  to  the  numerical  integration  step  size 
used  in  this  simulation  that  a  modification  is  necessary.  The  short 
correlation  time  and  hence  the  wideband  nature  of  this  random  variable 
may  easily  be  approximated  as  a  white  noise  input.  Thus,  the  user  clock 
error  model  used  in  this  study  is  given  in  Figure  12. 


u(t) 

Fig.  12.  Simplified  User  Clock  Error  Model  for  System 
The  suggested  values  for  CG,  Clt  and  C2  are  used;  the  initial  PMS 
value  of  the  covariance  on  the  white  noise  input  u(t)  is  assumed  to  be 
6.25  ft/sec. 

Thus , 
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and 

q22  -  6.25  ft2 /sec  (106) 

When  the  satellite  signal  passes  through  the  ionosphere,  the  signal 
is  bent  or  refracted.  Thus,  the  signal  path  is  not  a  perfectly  straight 
line  but  has  a  slight  amount  of  curvature  to  it. 


-O  —  Satellite 


indicated, 
range 


user 


True  Range 


Elevation  Angle 


Fig.  13.  Ionospheric  Delay 

This  increase  in  signal  path  length  due  to  unpredictable  ionospheric 
delay  is  on  the  order  of  15  to  25  feet  depending  upon  the  angle  of 
elevation  of  the  satellite  with  respect  to  the  user.  A  detailed  error 
source  model  for  this  ionospheric  delay  as  applied  to  a  set  of  syn¬ 
chronous  satellites  is  given  in  Reference  1.  However,  it  is  felt  that 
this  particular  model  is  not  applicableto  the  non-synchronous  satellite 
case.  As  mentioned  in  Chapter  III  the  errors  due  to  ionospheric  delay 
will  be  assumed  included  in  the  phase/range  satellite  clock  errors  for 
the  purpose  of  this  report.  An  area  for  future  study  would  be  the  syn¬ 
thesis  of  an  accurate  error  source  model  for  this  ionospheric  delay  in 
the  utilization  of  non-synchronous  satellites. 

The  final  error  sources  used  in  the  system  model  are  the  range 
biases.  These  are  modeled  simply  as  arbitrary  initial  conditions  on 
integrators  and  account  for  the  effect  of  tropospheric  delay  and  un¬ 
certainty  in  the  speed  of  light.  The  cumulative  magnitude  of  these  two 


Menaiii 
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error  sources  is  on  the  order  of  10  feet.  This  requires  four  additional 
state  variables  (one  for  each  satellite)  in  the  system  state  vector  but 
adds  no  additional  nonzero  terms  to  the  F  or  O'  matrices. 


System  Veasurerent  Equation 

The  measurement  equation  consists  of  a  set  of  four  range  divergence 
equations  as  derived  in  Chapter  111  (Equation  33) .  Three  of  these  are 
needed  to  reduce  the  position  error  and  the  fourth  is  required  to 
synchronize  the  user  clock.  This  study  is  concerned  with  us inf  only 
range  measurements  in  the  filter  design.  Incorporating  an  altimeter 
would  require  an  additional  measurement.  If  range-rate  measurements 
were  taker,  four  more  equations  would  be  required. 

The  system  measurement  equation  is 


h85b+^ 


(107) 


where  the  44  variables  of  x„  are  listed  in  Table  IV  and 

Ari 

Ar2 
*r3 


(108) 


is  the  4  by  44  matrix  shown  in  Figure  14.  The  white  Gaussian  noise 
vector  v  is  described  by  the  covariance  matrix  R  where 


Rs  " 


0 

0 

0 


0  0  0 
0 


o2  0 


0  a2  0 
r 


(109) 
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The  value  or  is  a  measure  of  the  white  noise  corrupting  the  high  fre¬ 
quency  satellite  signal  and  a  typical  value  of  or  -  10  feet  is  used  in 
this  design. 
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VI.  Filter  Designs 

This  chapter  discusses  a  few  of  the  various  Kalman  filter  designs 
which  were  simulated  using  the  44  state  reference  system  of  Chapter  V. 

Not  all  of  the  filter  designs  which  were  analyzed  are  presented;  only 
those  which  were  significant  and  illustrate  important  design  differences 
are  outlined.  These  Include:  (1)  the  44  state  optimal  filter,  (2)  a 
fully-coupled  15  state  filter,  (3)  a  15  state  filter  with  all  weak 
coupling  terms  removed  from  the  INS  error  model,  and  (4)  a  ten  state  fil¬ 
ter  illustrating  degraded  performance  caused  by  eliminating  too  many 
states  from  the  filter. 

Selection  of  the  "best"  Kalman  filter  design  Involves  a  trade-off 
study  between  desired  accuracy  and  the  computational  time  required  by 
an  onboard  computer.  Obviously,  the  most  accurate  filter  is  the  optimal 
44  state  filter;  however,  the  burden  this  would  place  on  the  computer  is 
unacceptable.  On  the  other  hand,  the  ten  state  filter  would  substantially 
reduce  the  computational  time  but  yields  unacceptable  accuracy.  There¬ 
fore,  the  procedure  of  this  study  was  to  successively  eliminate  a  few 
states  at  a  time  from  the  filter  state  vector  and  simulate  each  of  these 
reduced-order  or  suboptimal  filter  designs  on  the  computer.  So  long  as 
the  accuracy  of  the  filter  performance  was  degraded  only  slightly  in  each 
step,  this  process  was  continued  until  a  minimum  acceptable  filter  state 
vector  dimension  was  reached. 

The  Optimal  Filter 

The  optimal  filter  is  simply  an  exact  replica  of  the  system  model 

(i.e. ,  the  same  44  state  system  and  measurement  model  as  given  in 

I 

Chapter  V).  This  Kalman  filter  computes  the  optimal  gains  or  weighting 
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coefficients  from  the  filter  design  and  applies  them  to  the  system.  In 
the  optimal  filter,  the  values  of  the  state  variables  for  which  these 
weighting  coefficients  are  computed  are  identical  to  the  state  variable 
values  to  which  the  weighting  coefficients  are  applied  (the  system  model 
is  theoretically  assumed  to  be  exact).  This  optimal  44  state  filter  will 
yield  the  most  accurate  performance  possible  for  the  given  set  of  system 
state  and  measurement  equations. 

Plots  of  the  RMS  values  of  the  nine  plant  state  (position,  velocity, 
and  attitude)  error  covariances  for  the  optimal  filter  are  shown  in 
Figures  15  through  23.  These  plots  provide  the  basis  against  which  all 
sub-optimal  filter  performance  will  be  compared.  Note  that  the  time  of 
flight  is  one  hour  (3600  seconds)  for  these  and  all  plots  in  this  report. 
Also,  the  vertical  axes  are  automatically  scaled  in  the  computer  plotting 
routine  and  may  vary  from  one  case  to  the  next.  For  example,  the  optimal 
filter  z-RMS  velocity  error  is  scaled  from  0.14  ft/sec  to  0.78  ft/sec, 
whereas  the  plot  of  this  sane  variable  for  the  15  state  filter  presented 
later  in  this  chapter  is  scaled  from  0.32  ft/sec  to  0.96  ft/sec. 

Fully-Coupled  15  State  Filter 

The  44  state  optimal  filter  was  gradually  reduced  by  eliminating 
the  least  significant  variables  of  the  state  vector.  No  serious  degrada¬ 
tion  of  performance  was  noted  as  this  state  vector  was  reduced  to  the 
15  state  model.  However,  filter  doL-'^.ns  of  lower  dimension  did  diverge 
significantly,  one  example  of  this  1:  ^ iven  in  the  final  section  of 
this  chapter.  The  intermediate  desi  . ;  between  the  44  state  optimal 
and  the  15  state  sub-optimal  f i 1  ml  very  little  variation  in 

performance.  In  fact,  the  increased  accuracy  provided  by  retention  of 
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any  of  the  state  variables  in  addition  to  the  15  selected  is  so  slight 
that  it  is  felt  filter  designs  of  higher  dimension  are  unwarranted 
(except  for  the  addition  of  one  or  two  states  when  an  altimeter  is  in¬ 
corporated  in  the  system).  Filter  designs  whose  state  variable  vectors 
were  of  lower  dimension  than  this  15  state  design  exhibited  significantly 
degraded  accuracy  in  performance  and  are  not  recommended. 

All  sub-optimal  filter  designs  presented  in  this  report  contain  a 
model  of  the  INS  error  states  (position,  velocity,  attitude).  These 
states  are  required  to  propagate  the  prediction  of  the  predictor-correc¬ 
tor  behavior  of  the  filter.  These  INS  plant  error  states  comprise  the 
first  nine  variables  of  the  filter  state  vector  which  is  shown  in 
Table  VI. 

The  IiiS  error  source  models  (gyro  drift  and  accelerometer  errors) 
required  12  state  variables  in  the  system  model.  These  can  be  reduced 
to  six  state  variables  by  approximating  the  linear  combination  of  a 
random  bias  and  an  expotentlally  correlated  random  variable  as  a  random 
walk  (Ref  12:11-14  thru  11-16).  However,  all  of  the  above  12  state 
variables  can  be  eliminated  from  the  filter  design  by  approximating 
these  error  sources  as  white  Gaussian  noise  on  the  plant  states.  The 
accelerometer  error  sources  are  modeled  as  equivalent  driving  white 
noises  on  the  velocity  error  differential  equations  and  the  gyro  error 
sources  are  modeled  as  equivalent  white  noise  Inputs  on  the  attitude 
error  differential  equations. 

These  equivalent  white  noises  are  calculated  as  qj  values  which 
provide  approximately  the  same  amount  of  RMS  error  build-up  over  the 
given  interval  of  time  (3600  sec)  as  the  system  error  sources  which 
have  been  discarded  in  the  filter  design.  This  may  be  done  with  little 


66 


GA/EE/74M-1 


Table  VI 

Filter  15  State  Vector  Definition 

Error 


State 

Symbol 

Definition 

INS  Plant  Error  States 

1 

Ax 

x  position  error 

2 

Ay 

y  position  error 

3 

Az 

z  position  (altitude) 

4 

Ax 

x  velocity  error 

5 

Ay 

y  velocity  error 

6 

Az 

z  velocity  error 

7 

♦x 

x  attitude  error 

8 

<f>y 

y  attitude  error 

9 

tyz 

z  attitude  (heading)  error 

'Satellite  and  User  Clock  Errors 

10 

6TU 

user  clock  phase/range  error 

11 

X11 

user  clock  frequency  offset 

12 

6TSl 

sat.  clock  01  phase/range 
error 

13 

«TS2 

sat.  clock  02  phase/range 
error 

14 

6TS3 

sat.  clock  if 3  phase/range 
error 

15 

6Ts4 

sat.  clock  i/4  phase/range 

error 


RMS  Initial 
Condition 

3000  ft 
3000  ft 
300  ft 
2.0  ft/sec 
2.0  ft/sec 
0.5  ft/sec 
0.14  milllrad 
0.14  milllrad 
2.0  milllrad 

10,000  ft 
5  ft/sec 
100  ft 

100  ft 

100  ft 

100  ft 
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!  O 

i 

j 


loss  In  performance  so  long  as  the  measurement  sample  frequency  is  high  com¬ 
pared  with  dominant  system  time  constants  (Ref  1:4-3  thru  4-4).  This 
requires  six  diagonal  elements  in  the  filter  Q'  matrix  corresponding  to 
these  equivalent  white  noise  inputs.  The  general  equation  used  to  cal¬ 
culate  these  q*  values  is 


qj 


* k. 

Si' 


20biTi 


(110) 


where  is  the  initial  covariance  of  the  corresponding  system  bias 
error  source  and  Si  and  Ti  are  the  inverse  correlation  time  and  correla¬ 
tion  time,  respectively,  of  the  corresponding  system  expotentlally  corre¬ 
lated  error  source.  The  three  equivalent  noises  on  the  velocity  error 
differential  equations  are 


q*  ■  2o2  t 
*»  abx  arx 

q*  -  (2) [(66  x  10"6  g) (32.2  ft/sec2-g) ]2 (600  sec) 
q*  ■  5.411  x  10"3  ft2/sec3 


(111) 

(112) 

(113) 


Similarly, 

q$  ■  q$  (114) 

q*  ■  2.164  x  10”2  ft2/sec3  (115) 

0 

The  three  equivalent  noises  on  the  attitude  error  differential  equations 
are 


<1? 


2o2  T 
ebx  erx 


(116) 


q*  "  (2)  [(0.025  deg/hr)  (4. 848  x  10-6  )  ]2  (3600  sec)  (117) 


O 


q*  ■  1.0576  x  10“3®  rad2 /sec 


(118) 
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Similarly, 

q*  -  q*  (119) 

q*  -  2.7076  x  10“10  rad2/sec  (120) 

iiote  that  since  there  are  no  white  noise  inputs  on  the  position  differ¬ 
ential  equations 

q*  -  q*  -  q*  -  0  (121) 


The  four  range  bias  error  sources  (on  the  order  of  10  ft)  are  easily 
eliminated  by  simply  increasing  the  value  of  the  diagonal  elements  of 
the  measurement  noise  matrix  (R)  to  20  feet. 

Upon  examination  of  the  clock  error  source  models  and  the  relative 
magnitudes  of  their  state  variable  initial  conditions  as  presented  in 
Chapter  V,  it  is  apparent  that  the  most  significant  variables  are  the 
four  phase/range  errors  of  the  satellite  clocks  (in  the  order  of  100  ft) , 
the  user  clock  phasa/range  error  (on  the  order  of  10,000  ft),  and  the 
user  clock  frequency  offset  error  (on  the  order  of  5  ft/sec).  These 
six  variables  are,  in  fact,  the  only  ones  used  in  this  15  state  filter 
design  except  for  the  nine  IUS  plant  error  states. 

The  user  clock  is  simplified  by  dropping  the  long  term  stability 
state.  The  filter  model  for  the  user  clock  is  shown  in  Figure  24. 
Initial  conditions  and  noise  values  remain  unchanged  from  the  system 
model. 

$TU 

Fig.  24.  Simplified  User  Clock  Error  Model 
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As  in  the  system,  the  noise  value  input  to  the  phase  integrator  is 

q*0  -  6.25  ft2 /sec  (122) 

The  error  source  model  for  the  four  satellite  clocks  is  greatly 
simplified  by  eliminating  all  satellite  clock  variables  from  the  filter 

state  vector  except  for  the  phase/range  errors.  The  simplified  sate¬ 
llite  clock  error  model  for  the  15  state  filter  is  given  in  the  follow¬ 
ing  figure. 


Fig.  25.  Simplified  Satellite  Clock  Error  Model 
Although  the  satellite  clock  frequency  offset  state  has  been  de¬ 
leted  from  the  filter  state  vector,  the  phase  error  build-up  caused  by 
this  offset  cannot  be  ignored.  This  effect  is  accounted  for  by  finding 
the  white  noise  input  to  the  clock  phase  integrator  which  causes  an 
equivalent  build-up  in  phase  error  over  a  specified  time.  From  Equation 
(96)  it  is  seen  that  a  frequency  offset  of  Cj  will  cause  an  RMS  phase 
error  build-up  of 

6Tg  -  CjAt  (123) 

Similarly,  a  white  noise  input  of  qAt  to  the  phase  integrator  will  cause 
an  RHS  build-up  of 

6TS  -  (  t)1/2  (124) 

Substitution  of  the  right-hand  ;  '  3tion  (123)  into  the  left- 

hand  side  of  Equation  (124)  and  solv-ui^  yields  an  equation  for  the 
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required  equivalent  q  value  of  the  Q1  matrix  (Ref  1:4-15  &  4-16). 

q  -  C2  At  (125) 

Tills  simple  satellite  clock  model  need  be  valid  only  over  the 
satellite  clock  update  interval,  or  flight  time.  Thus  At  was  set  at 
one  hour  (3600  sec). 

q*2  -  (.01)2(3600)  -  0.36  ft2/sec  (126) 

Similarly, 

’*3  "  ■  “Is  "  ‘Jta  <l27> 

The  F  matrix  representing  this  15  state  filter  design  is  given  in 
Equation  (128).  The  "Pinson  Error  Model"  is  the  same  as  given  in 
Figure  7. 


1-9 

10-11 

12-15 

Pinson  1 

1 

— 

1-9 

Error  . 
Model 

0  I 

~o7  '■ 

0 

10-11 

0  | 

i 

-°-i 

0  I 

0 

12-15 

0  - 

0 

The  Q' ,  H,  and  R  matrices  are  identical  to  those  of  the  filter  design 
presented  next,  and  are  not  displayed  here  for  the  sake  of  brevity. 
Also,  the  RMS  error  plots  for  this  filter  design  were  indistinguishable 
from  those  of  the  decoupled  version  discussed  next  and  are  omitted. 

Tills  filter  design  and  all  others  in  this  chapter  were  simulated 
using  30  second  measurement  rate.  In  the  subsequent  chapter,  this  15 
state  filter  is  used  as  a  basis  of  comparison  for  various  measurement 
rates. 
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15  State  Filter  -  Weak  Coupling  Removed 

This  filter  provided  the  best  tradeoff  between  accuracy  and  compu¬ 
tational  time  for  the  particular  problem  of  this  study  and  is  therefore 
the  recommended  Kalman  filter  design. 

Incorporating  the  Pinson  INS  error  model  into  the  filter  design  re¬ 
quires  the  computation  of  36  non-zero  terms,  in  the  F  matrix  many  of 
which  are  complex  (see  Figure  7).  However,  for  an  aircraft  cruising  at 
subsonic  speeds,  as  in  the  flight  profile  of  this  study,  many  of  these 
terms  are  negligible  and  may  be  deleted  from  the  F  matrix.  This  includes 
terms  on  the  order  of  magnitude  of  wie  such  as  Vx/R  and  Vy/R.  However, 
the  Schuler  frequency  terms  (u£  and  2u£)  are  included  so  that  the  model 
applicability  will  not  be  restricted  to  short  periods  of  time  (on  the 
order  of  30  minutes).  The  simplified  F  matrix  along  with  the  Q1,  H,  and 
R  matrices  for  this  15  state  filter  design  with  weak  coupling  terms  re¬ 
moved  are  given  in  Equations  (129)  through  (132). 
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The  values  for  q4  through  qJ5  are  as  given  in  the  preceding  section  of 
this  chapter. 


K 


f 


1234-9 


1 

2 

3 

4 


Arxl  iryl  irzl 
irx2  *ry2  irz2 
*rx3  *ry3  *rz3 
*rx4  iry4  irz4 


10 

1 

1 

1 

1 


11  12-15 

0  I  -1  0  0  0 

I 

o  0  -1  0  0 


0  |  o  0  -1  0 

0  *0  0  0  -1 

I  _ 


(131) 


Where  of  ■  20  ft. 

The  system  RMS  position,  velocity,  and  attitude  error  plots  for 
this  filter  design  are  given  in  Figures  26  through  34  and  should  be  com> 
pared  with  those  for  the  optimal  filter.  Table  VII  following  these 
plots  compares  the  system  plant  and  user  clock  phase/range  error  values 
for  both  the  optimal  and  the  15  state  decoupled  filter  at  2400  seconds 
with  the  initial  values  to  provide  a  quick  comparison  of  performance. 
These  are  "average"  values  as  explained  in  the  next  chapter. 

Note  that  the  effectiveness  of  both  the  optimal  and  the  15  state 
sub-optimal  filter  is  approximately  the  same  in  controlling  the  x,  y,  z 
position  and  user  clock  phase/range  errors.  Also,  the  15  state  filter 
yields  attitude  accuracy  close  to  that  of  the  optimal  filter.  However, 
the  x  and  y  velocity  errors  are  reduced  to  approximately  10%  of  their 
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Table  VII 


Average  RMS  Errors  of  INS  Plant:  States 
15  state  filter  and  optimal  44  state  filter  performance  compared 
against  Initial  values.  User  clock  phase/range  Is  also  included. 


Variable 

Initial  Condition 

Optimal  Filter 
(t  ■  2A00  sec) 

15  State  Filte: 
(t  ■  2A00  sec) 

Ax 

3000  ft 

A5.85  ft 

A7.A7  ft 

Ay 

3000  ft 

125.2  ft 

126.1  ft 

As 

300  ft 

256.3  ft 

264.6  ft 

Ax 

2  ft/aec 

0.103  ft/8ec 

0.213  ft/sec 

Ay 

2  ft/sec 

0.117  ft/8ec 

0.217  ft/sec 

As 

0.5  ft/sec 

0.158  ft/sec 

O.AAA  ft/sec 

♦x 

0.1A  milllrad 

0.065  milllrad 

0.082  milllrad 

♦y 

0.1A  milllrad 

0.065  milllrad 

0.08A  milllrad 

♦* 

2  milllrad 

1.27  milllrad 

1.A5  milllrad 

«TU 

10,000  ft 

199.0  ft 

205.6  ft 

83 


GA/EE/74M-1 


initial  values  by  the  sub-optimal  filter  as  compared  to  about  5Z  for  the 
optimal  case.  The  vertical  velocity  error  is  reduced  only  slightly  be¬ 
low  the  initial  value  in  the  sub-optimal  case,  whereas  it  is  reduced  to 
approximately  1/3  its  initial  value  by  the  optimal  filter.  Also  note 
that  the  altitude  error  is  merely  kept  from  growing  above  its  initial 
value  in  both  cases  and  altitude  accuracy  is  not  significantly  improved. 
These  results  suggest  two  areas  for  future  investigation.  First,  range- 
rate  measurements  could  be  effective  in  improving  upon  velocity  errors; 
and  secondly,  the  addition  of  altimeter  measurements  would  almost  cer¬ 
tainly  reduce  the  magnitude  of  the  altitude  error.  However,  x  and  y 
position  error,  which  is  usually  of  primary  importance,  is  very  effective¬ 
ly  controlled  by  this  13  state  filter  design. 

10  State  Filter 

A  filter  design  utilizing  only  ten  states  is  presented  here  for  the 
purpose  of  illustrating  degraded  performance  characterized  by  reducing 
the  dimension  of  the  filter  state  vector  too  far.  In  this  case,  the 

t 

Kalman  gains  computed  by  the  filter  do  not  include  the  effect  of  signi¬ 
ficant  variables  which  have  been  deleted  from  the  filter.  Thus,  when 
these  gains  are  applied  to  the  system,  the  divergent  growth  of  the  error 
equations  is  not  effectively  controlled. 

In  this  case,  the  nine  plant  error  state  variables  and  the  user 
clock  phase/range  error  are  the  only  ones  retained  in  the  filter.  Thus, 
this  design  is  obtained  from  the  13  state  filter  by  discarding  the  four 
satellite  clock  phase/range  error  states  and  the  user  clock  frequency 
offset  error.  The  decoupled  version  of  the  INS  error  model  (Equation 
(129))  is  used  in  the  F  matrix  and  the  user  clock  is  modeled  in  a  manner 
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similar  to  the  simplified  satellite  clock  model  of  the  15  state  filter. 
The  H  matrix  is  obtained  by  deleting  the  four  satellite  clock  terms 
(-l’s)  from  the  15  state  design.  Also,  the  RMS  value  of  the  measurement 
noise  (or)  is  increased  to  100  feet  to  account  for  the  increased  un¬ 
certainty  in  the  range  measurement  value.  The  RMS  error  plots  for  this 
ten  state  filter  design  illustrate  its  unacceptable  performance  and  are 
shown  in  Appendix  A. 
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VII.  Measurement  Update  Rate  and  Selection  of  Observables 

The  first  section  of  this  chapter  investigates  the  effect  of 
varying  the  measurement  update  rate.  The  fully  coupled  15  state  filter 
design  of  the  previous  chapter  is  used  in  this  investigation.  The  mea¬ 
surement  update  rate  requirements  may  vary  greatly  from  system  to  system 
depending  on  the  accuracy  required  as  well  as  the  capabilities  of  the 
equipment  used  to  take  the  measurements  for  the  update.  For  example, 
a  smaller  onboard  computer  would  normally  require  a  larger  computation¬ 
al  time.  Additionally,  a  comparison  is  made  between  obtaining  measure¬ 
ments  from  satellites  which  are  randomly  selected  from  those  in  view 
and  obtaining  measurements  from  a  sequence  of  "optimum"  sets  of  four 
satellites,  three  of  which  form  a  "most  nearly  orthogonal"  set  at  the 
time  of  measurement. 

The  update  rates  used  in  this  report  are  as  follows:  5  seconds, 

15  seconds,  30  seconds,  60  seconds,  and  90  seconds.  The  results  of 
these  update  rates  are  compared,  at  a  specified  time  in  the  flight 
profile,  in  Table  VIII  at  the  end  of  this  chapter. 

Comparison  of  Update  Rates 

A  few  representative  plots  of  the  RMS  errors  for  the  nine  INS 
plant  error  states  are  shown  on  the  following  pages.  In  the  previous 
chapter,  a  30  second  update  rate  was  used  in  determining  the  Kalman 
filter  state  vector  dimension.  ThcrJi'ore,  plots  of  the  nine  plant  error 
states  with  this  update  rate  are  Inc  jded  in  that  chapter  and  will  not 
be  repeated  here.  Plots  of  the  thro.'  error  state  variables;  x  position, 
y  attitude  and  z  velocity  for  uprf.-t  of  5  seconds,  15  seconds, 

60  seconds,  and  90  seconds  appear  in  i’igures  35  through  46.  Plots  of 
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the  remaining  six  plant  error  states  for  these  measurement  rates  are 
included  in  Appendix  B. 

It  should  be  noted  that  the  horizontal  axis  or  time  scale  remains 
constant  for  all  the  plots,  but  the  vertical  axis,  RMS  error,  was 
scaled  for  each  plot  to  effectively  show  the  results  of  the  update  rate. 
The  plotting  program  was  set  to  plot  points  at  a  minimum  interval  of 
30  seconds.  Because  of  this  plot  control  the  initial  updates  occurlng 
in  the  0  to  30  second  time  period  are  not  shown  for  the  5  and  15  second 
measurement  rates.  In  these  two  cases  the  errors  have  been  greatly  re¬ 
duced  by  the  time  the  first  points  are  plotted  and  the  lower  region  of 
the  vertical  axis  has  been  greatly  magnified.  This  illustrates  the 
divergent  behavior  of  the  position  error  states  at  the  end  of  this 
flight  profile.  This  aspect  will  be  discussed  in  the  final  section  of 
this  chapter. 

These  plots  indicate  that  increasing  the  update  rate  does  not 
greatly  improve  the  RMS  errors.  In  general,  a  slower  update  rate  allows 
the  errors  to  grow  to  higher  values  during  the  initial  portion  of  the 
flight  (first  15  minutes).  The  error  states  after  this  initial  time 
period  exhibit  an  almost  "steady  state"  time  behavior.  The  slower  up¬ 
date  rates  provide  R11S  error  values  which  are  a  little  higher  on  the 
average  than  the  faster  measurement  rates.  This  fact  is  illustrated  in 
Table  VIII.  This  table  shows  the  average  RMS  errors  of  the  plant  error 
states  at  each  update  rate,  for  a  pa  ticular  time  during  the  flight 
profile.  Also  shown  in  the  table  ir  m  per  cent  improvement  from  the 
90  second  update  rate  to  the  5  secon  update  rate.  "Average"  values 
are  obtained  by  fitting  a  smooth  i:n  ■  o  the  data  points.  This  is 
sccompllshed  by  adding  the  RMS  error  values  of  a  variable  at  the  time 
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instant  before  a  measurement  to  its  value  just  after  the  measurement 
and  dividing  by  two.  Also,  slower  update  rates  permit  the  errors  to 
grow  to  substantially  higher  values  between  measurements. 

Selection  of  Observables 

When  using  SAMUS  (a  computer  program  detailed  in  Chapter  VIII) ,  two 
test  cases  were  executed  using  a  twenty-two  state  optimal  filter.  In 
both  cases  four  satellites,  or  observables,  were  selected  for  each  re¬ 
set  segment  by  studying  data  output  of  the  profile  generator  program  of 
Chapter  III,  which  calculates  unit  vectors  in  the  navigation  frame  from 
the  user  to  each  satellite.  Three  of  these  range  measurements  sre  re¬ 
quired  for  correcting  x,  y,  and  z  position  and  the  fourth  measurement 
is  needed  for  synchronization  of  the  user  clock.  Three  satellites  hav¬ 
ing  the  largest  components  along  the  x,  y,  and  z  axes,  respectively, 
when  compared  with  all  other  satellites  ln-vlew  at  a  particular  time, 
are  chosen.  These  form  a  most  nearly  orthogonal 'set  of  measurements. 

The  fourth  satellite  is  chosen  as  the  one  having  the  second  largest  com¬ 
ponent  along  the  z  axis.  For  example,  the  "optimal"  set  of  satellites 
at  the  particular  time  instant  shown  in  Table  III  is  composed  of  sate¬ 
llites  52,  13,  93,  and  62. 

In  case  one,  the  flight  time  was  divided  equally  into  two  reset 
segments.  Since  the  satellite  unit  vectors  are  time  varying,  the  sat¬ 
ellites  chosen  for  each  reset  segment  provide  the  best  information,  on 
the  average,  over  the  given  time  interval.  SAMUS  was  then  executed, 
using  these  two  reset  segments,  to  calculate  the  values  of  the  IKS 
error  states  as  a  function  of  time  over  the  one  hour  flight  profile. 

The  data  from  this  run  was  next  input  to  a  general  covariance  plotting 
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program  that  provided  the  plots  for  the  root  mean  square  (RMS)  error  of 
the  diagonal  elements  of  the  covariance  matrix  for  the  twenty-two  error 
states.  Figure  47  is  a  plot  of  the  x  and  y  position  error ,  in  nautical 
miles,  for  the  time  period  of  the  profile.  Figure  48  is  a  plot  of  the 
altitude  error,  in  feet,  for  the  same  period.  The  results  for  the  posi¬ 
tion  error  indicate  a  decrease  from  .5  nautical  mile  initially  to  .005 
nautical  mile  after  the  first  reset.  The  error  remains  constant  for 
about  30  minutes  then  the  error  starts  to  Increase.  The  altitude  error 
indicate*  a  decrease  from  100  feet  initially  to  15  feet  after  the  first 
reset.  At  this  time  the  error  starts  to  increase  until  it  reaches  approx¬ 
imately  50  feet  at  the  end  of  the  flight  profile. 

In  case  two  of  the  same  22  state  optimal  filter  was  used;  however, 
this  run  contained  six  reset  segments.  Since  the  satellite  unit  vectors 
were  time  varying,  more  accurate  satellite  information  could  be  acquired 
by  changing  the  observables,  satellites,  at  shorter  intervals.  Figure 
49  shows  the  plot  of. the  x  and  y  position  error  for  this  run.  The  re¬ 
sults  show  that  the  position  errors  did  not  diverge  at  the  end  of  the 
profile.  Figure  50  Illustrates  the  elimination  of  divergent  behavior 
in  the  altitude  error. 

A  comparison  of  the  two  cases  indicates  that  by  increasing  the  num¬ 
ber  of  reset  segments  so  that  measurements  are  taken  from  "optimal"  sets 
of  satellites  throughout  the  flight,  the  IHS  plant  errors  reach  almost 
steady-state  values  and  the  tendency  toward  divergent  behavior  in  the 
latter  part  of  the  flight  disappears.  However,  the  general  covariance 
program  used  in  generating  the  bulk  of  the  plots  appearing  in  this  re¬ 
port  does  not  have  the  capability  for  multiple  reset  segments.  Thus,  a 
aet  of  4  satellites  giving  reasonably  good  performance  throughout  the 
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one  hour  flight  was  used  in  the  simulation  work  of  Chapters  VI  and  VII. 
Therefore,  the  divergent  behavior  illustrated  in  Figures  35  and  36  can 
probably  be  eliminated  by  obtaining  measurements  from  a  sequence  of 
"optimal"  sets  of  satellites.  Incorporating  this  capability  into  the 
present  covariance  analysis  program  is  an  area  suggested  for  future  study. 
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VIII,  Computer  Simulation  Programs 

This  chapter  will  outline  the  computer  programs  utilised  in  the 
simulation  of  an  INI  flight.  The  first  section  is  a  discussion  of  the 
State  space  Analysis  of  MUltisensor  Systems  (SAMUS)  program.  SAMUS  was 
used  in  the  early  stages  of  this  design  study,  and  later  abandoned  in 
favor  of  a  set  of  two  computer  programs  which  will  be  discussed  in  the 
second  section.  This  second  set  of  computer  programs,  consisting  of  a 
flight  profile  generator  and  a  covariance  analysis  program,  forms  the 
analysis  tool  used  in  this  study. 

SAMUS 

The  first  computer  program  utilized  was  the  State  space  Analysis  of 
MUltisensor  Systems  (SAMUS).  This  is  a  general  purpose  computer  program 
for  performing  error  analysis  of  cruise  Inertial  navigation  systems. 

SAMLiS  was  developed  and  documented  (kef  6)  under  a  computer  aided  design 
effort  funded  by  North  American  Rockwell  Corporation. 

The  program  was  written  to  allow  the  user  to  perform  analytical  statis¬ 
tical  error  analysis  of  typical  cruise  system  problems  with  little  or  no 
programming  required.  It  has  been  preprogrammed  to  give  the  equations  of 
motion,  a  master  system  description  matrix,  a  master  measurement  matrix, 
aud  the  error  source  statistics  and  coordinate  frames  which  are  normally 
used.  In  a  typical  cruise  system  utilizing  a  Kalman  filter,  estimation 
and  control  occur  at  short,  regular  i,  tervals.  This  precludes  specify¬ 
ing  each  control  occurrence  Individ*'.  -  Thus,  the  program  was  written 

so  that  the  resets  could  be  specified  sepments.  Each  observable  sneci- 

fied  in  a  segment  has  an  initial  c  :  n  tine  in  the  segment  and  Is 

observed  at  a  regular  Interval  for  the  rest  of  the  reset  segment. 
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Cruise  system  problems  require  investigation  of  a  variety  of  error 
source  statistics.  The  SAMUS  program  allows  the  user  to  specify  error 
statistics  ranging  from  a  constant  to  a  second  order  noise.  The  program 
has  the  system  distribution  matrices  of  the  typical  cruise  error  sources 
preprogrammed.  The  user  can  specify  required  error  sources  by  means  of 
an  error  source  code  list  (Ref  6:1-2). 

however,  the  program  was  not  designed  to  provide  the  necessary 
measurement  error  sources  that  are  needed  for  the  problem  undertaken  in 
this  report.  In  the  system  under  investigation  the  measurements  are 
taken  from  satellites.  A  model  of  a  clock  was  necessary  for  the  user 
and  each  satellite.  A  major  modification  of  the  program  would  have  been 
required  to  model  the  error  sources  for  these  clocks.  Also,  the  computa¬ 
tional  burden  of  5AMJS  was  prohibitive;  for  example  it  required  approx¬ 
imately  one  hundred  and  forty  thousand  octals  of  computer  memory  and 
six  hundred  seconds  of  computer  time  to  perform  the  calculations  for  a 
simplified  22  state  filter  pass.  It  was  for  these  reasons  that  a  more 
adequate  and  efficient  program  was  needed  for  this  Investigation. 

Computer  Program  Final  Selection 

A  combination  of  two  separate  computer  routines  was  adopted  to 
perform  the  required  analysis  of  this  report.  The  first  program  is  a 
flight  profile  generator.  This  program,  given  the  necessary  initial 
conditions,  computes  the  parameters  of  the  flight  profile.  This  infor¬ 
mation  is  then  stored  in  a  permanent  file.  The  second  program  is  a  gener¬ 
al  covariance  program.  This  program  propagates  the  system  and  filter 
state  variable  covariances  forward  in  time  by  numerical  integration  tech¬ 
niques,  reading  from  the  flight  profile  the  values  for  aircraft  latitude. 
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longitude,  altitude,  range,  velocity,  acceleration  and  satellite  posi¬ 
tion  information.  These  values  are  calculated  in  the  flight  profile 
generator  at  the  time  Instant  corresponding  to  the  numerical  integra¬ 
tion  step  size  of  the  general  covariance  program.  The  filter  covari¬ 
ances  are  propagated  forvard  until  the  update  time  is  reached.  The 
system  is  then  propagated  forvard  to  the  update  time.  Because  the  fil¬ 
ter  and  the  system  are  propagated  in  a  parallel  manner  requiring  dupli¬ 
cate  information,  two  identical  flight  profiles  vere  generated.  Logic 
was  incorporated  in  the  general  covariance  program  to  read  data  from 
the  first  flight  profile  for  the  filter  propagation  and  from  the  second 
flight  profile  for  the  system  propagation.  By  precalculating  the  flight 
profile  and  storing  the  data,  computer  memory  and  time  vere  saved. 

Profile  Generator.  The  flight  profile  generator  (PROFGEN)  was 
written  to  compute  the  position,  velocity  and  acceleration  of  a  point 
mass  that  is  moving  above  the  earth's  surface.  The  user  must  supply  the 
program  initial  values  for  position,  attitude,  and  velocity  plus  commands 
to  turn,  fly  straight,  change  heading  in  sinusoidal  fashion,  accelerate 
or  decelerate.  The  program  puts  out  position  data  including  latitude, 
longitude,  altitude  and  range  from  the  earth's  center  and  velocity  and 
acceleration  data  in  the  North-West-Up  navigation  coordinate  frame. 
Additionally,  the  flight  profile  generator  was  modified  to  calculate  the 
latitude,  longitude,  In-vlev  criterion,  and  unit  vector  components  for 
each  of  the  27  satellites  at  every  specified  integration  time  step. 

This  portion  of  the  flight  profile  generator  comprises  the  satellite 
motion  generator  discussed  in  Chapter  III. 

The  investigation  undertaken  in  this  report  used  one  of  the  two 
straight  flight  options  available  in  the  program.  This  option  was  a 
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sector  of  a  great  circle  path  around  the  earth.  The  great  circle  path 
varies  the  heading  angle  to  maintain  flight  In  a  fixed  plane  passing 
through  the  earth's  center.  This  option  was  used  with  a  constant  veloc¬ 
ity  for  a  straight  and  level  flight  of  one  hour  duration.  Table  IX 
specifies  the  values  for  the  flight  profile  used  in  this  investigation. 

Table  IX 
Flight  Profile 


Initial  Head 

45* 

Initial  Latitude 

30*N 

Initial  Longitude 

75  *W 

Altitude 

10,000  ft. 

Velocity 

600  mph 

Duration  of  Flight 

60  minutes 

Note  that  the  above  parameters  must  first  be  projected  onto  the  naviga¬ 
tion  frame  (for  example ,  Vn  -  [Vx,  Vy,  VZ]T)  before  they  are  utilized  in 
the  proolem  solution. 

The  input  data  that  the  user  supplies  to  execute  this  program  speci¬ 
fies  the  initial  conditions,  the  flight  maneuvers  of  the  aircraft,  and 
the  method  of  solution  of  the  differential  equations.  The  program  uses 
a  NAMELIST  data  input  format  that  permits  the  entry  of  character  strings 
consisting  of  parameter  names  with  their  values  in  the  user's  choice  of 
format  specification.  Two  NAMELIST  :  ut  lists  are  used,  PRDATA  and 
PASDATA.  The  PRDATA  (Problem  Data)  1  :1s  called  at  the  start  of  each 
problem.  This  list  contains  nine  pro.  ism  parameters,  all  of  which  must 
be  specified  and  none  of  which  ch..  .  Lag  execution  of  the  problem. 

Ilia  purpose  of  these  parameters  is  to  specify  all  initial  conditions. 
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Values  from  the  PASDATA  (Path  Segment  Data)  list  are  called  at  the  start 
of  each  path  segment.  This  contains  twelve  path  parameters  that  specify 
the  maneuver  to  be  followed  in  a  segment  and  the  method  of  integration 
to  use.  A  problem  may  be  divided  into  a  maximum  of  20  path  segments 
(Ref  13:1-5). 

General  Covariance  Program.  The  second  program  used  in  this  investi¬ 
gation  was  a  general  covariance  program.  This  program  propagates  the 
system  and  filter  state  variables  forward  in  time  by  numerical  integra¬ 
tion  techniques.  There  are  four  subroutines  in  this  program  that  the 
user  must  rodify  for  the  particular  system  under  study.  The  first  sub¬ 
routine  to  be  modified  (FLMAT)  supplies  the  filter  matrix  elements.  Each 
non  zero  element  of  the  F  matrix  is  input  as  either  a  constant  or  a 
function  of  time.  A  set  of  indices  for  each  non-zero  element  is  then 
specified.  The  indices  are  used  to  locate  the  numbered  elements  in  the 
F|  matrix,  for  example;  the  element  A(4)  may  have  Indices  3,4  which 
would  locate  it  in  the  third  row  and  the  fourth  column  of  the  Ff  matrix. 
The  Hf  matrix  was  specified  the  same  way. 

The  system  matrix  (SYSMAT)  is  the  second  subroutine  to  be  modified. 
This  routine  is  set  up  the  same  as  TLtiAT,  specifying  only  the  non-zero 
elemeuts  of  the  Fa  and  Hs  matrices. 

The  next  subroutine  to  be  changed  is  Trajectory  (TRAJ).  This  sub¬ 
routine  is  called  at  each  Integration  time  Instant  and  reads  the  flight 
profile  Information  from  one  of  the  two  permanent  files,  depending  on 
the  type  of  pass  being  made,  filter  or  system.  The  profile  generator 
provides  the  values  for  aircraft  latitude,  longitude,  altitude,  range, 
velocity,  acceleration  and  satellite  position  information.  The  angular 
rates  and  angular  accelerations  of  the  coordinate  frames  which  are  used 
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to  form  the  Pinson  error  model  are  then  calculated  using  this  informa¬ 
tion.  This  data  was  placed  in  a  labeled  common  statement  so  that  it 
could  be  used  in  the  other  subroutines  requiring  these  values. 

The  last  subroutine  to  be  modified  is  User  Input  (USRIN)*  The 
values  for  Q  and  K  are  specified  for  both  the  filter  and  the  system. 

This  information  was  also  placed  in  a  labeled  common  statement  so  that 
the  values  could  be  used  in  FLMAT  and  SYSMAT  subroutines. 

This  program  proved  to  be  much  more  efficient  than  SAMUS.  The  22 
state  filter  that  was  executed  using  SAHUS  required  approximately  one 
hundred  and  fourty  thousand  octals  of  computer  memory  and  six  hundred 
seconds  of  computer  time.  The  general  covariance  program  required  only 
forty-five  thousand  octals  of  memory  and  tvo  hundred  seconds  of  computer 
time  to  execute  the  same  22  state  filter  pass.  The  efficiency  of  this 
program  comes  primarly  from  the  fact  that  it  only  performs  operations 
on  the  non  zero  elements  of  the  F,  H,  Q,  and  R  matrices.  This  saves 
computer  memory  as  wqll  as  execution  time.  For  example,  the  forty-four 
state  system  model  would  require  nineteen  hundred  octals  of  memory  for 
the  entire  matrix.  However,  there  were  only  seventy-three  non  zero 
elements  in  this  matrix.  This  resulted  in  a  savings  of  over  eighteen 
hundred  octals  of  memory. 

Covariance  Plotting  Program.  Once  the  general  covariance  program 
has  calculated  the  diagonal  elements  of  the  covariance  for  each  state 
at  a  particular  time  Instant,  they  are  read  into  a  permanent  file.  The 
plotting  program  reads  these  values  from  this  file  and  prepares  them  for 
off-line  plotting.  This  information  is  then  placed  on  a  magnetic  tape 
and  can  be  plotted  at  any  time.  This  program  plots  the  RMS  error  versus 

I 

time  for  the  system  and  the  filter  over  the  sixty  minute  flight  profile. 


113 


GA/EE/74M-1 


The  program  is  capable  of  plotting  the  RMS  v clues  on  a  regular  linear 
scale  or  a  log  scale. 

A  brief  description  of  the  computer  tools  used  in  this  investiga¬ 
tion  was  presented  in  this  chapter.  It  was  shown  that  the  computation¬ 
al  burden  of  the  original  computer  program  (SAMUS)  was  prohibitive. 

The  programs  used  in  the  final  analysis  proved  to  he  much  more  efficient. 
A  listing  of  these  programs  was  not  Included  in  this  report  because  of 
their  extensive  length. 
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IX.  Conclusions  and  Recommendations 


Conclusions 

Based  on  tne  material  presented  in  this  thesis,  as  well  as  knowledge 
gained  throughout  this  study,  the  following  conclusions  are  drawn: 

1.  The  15  state  filter  vector  consisting  of  nine  INS  error  states, 
two  user  clock  errors,  and  four  satellite  clock  errors,  provides  the 
best  performance  tradeoff  between  high  accuracy  and  low  computational 
time. 

2.  The  15  state  filter  with  weak  coupling  terms  removed  from  the 
Pinson  IKS  error  model  provides  almost  identical  performance  as  the 
fully-coupled  filter  for  this  particular  flight  profile.  This  simplifi¬ 
cation  results  in  a  substantial  reduction  of  the  computational  burden  on 
the  user's  computer. 

3.  Measurement  rates  of  once  every  5,  15,  30,  60,  or  90  seconds  all 
provide  satisfactory  filter  performance  (depending  upon  the  user's  re¬ 
quirements).  The  slower  update  rates  provide  slightly  larger  RMS  errors, 
on  tne  average,  and  allow  more  error  build-up  between  measurements. 

4.  Frequent  changing  of  satellite  observables  so  that  measure¬ 
ments  are  taken  from  "optimal"  sets  of  satellites  (as  explained  in 
Chapter  VII)  prevents  divergence  of  the  error  states  in  the  latter  part 
of  the  flight. 

5.  The  vertical  velocity  and  altitude  errors  are  the  most  diffi¬ 
cult  to  control.  Use  of  an  altimeter  in  conjunction  with  the  filter 
should  eliminate  this  divergent  tendency. 

6.  Acceptable  accuracy  of  the  INS  errors  is  obtained  in  this 
study  using  only  range  measurements.  The  simulation  results  suggest 
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that  Inclusion  of  range-rate  measurements  could  possibly  increase  the 
filter  accuracy  (especially  in  the  velocity  errors)* 

7.  Finally*  it  should  again  be  emphasized  that  this  type  of  co- 
variance  analysis  is  very  sensitive  to  (a)  the  arbitrary  initial  con¬ 
ditions  of  the  covariance  matrix  P(0)  (in  the  initial  transient  por¬ 
tion),  (b)  erroneous  system  reference  models,  and  (c)  the  dynamics  of 
the  flight  profile. 

Recommendations 

The  following  recommendations  are  made  for  continued  study  in  the 
subject  area  of  this  thesis  and  for  Improvement  of  the  existing  computer 
program: 

1.  Incorporation  of  an  altimeter  into  the  filter  and  system  should 
provide  more  effective  control  of  the  vertical  velocity  and  altitude 
errors*  This  would  require  the  addition  of  two  or  three  INS  plant  error 
states  in  both  the  reference  system  and  the  filter  state  variable  vector. 
Additionally,  this  would  add  one  equation  to  the  syaten  and  filter  set 
of  measurement  equations. 

2*  Inclusion  of  an  accurate  linear  model  of  the  ionospheric  de¬ 
lays  would  Improve  the  reference  system's  accuracy  in  describing  the 
"real  world"  problem.  A  model  applicable  to  the  non-synchronous  sate¬ 
llite  case  night  be  derived  through  modification  of  the  ionospheric  de¬ 
lay  model  presented  in  Reference  1  for  the  synchronous  case. 

3.  The  effect  upon  filter  performance  of  including  range-rate 
measurements  should  be  Investigated.  A  comparison  with  the  usage  of 
only  range  measurements  could  then  be  drawn.  Use  of  satellite  bearing 
measurements  might  also  be  studied. 
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4.  An  optimal  sequencing  scheme  should  be  developed  to  automatically 
select  the  best  satellites  from  which  to  take  measurements  of  those  in¬ 
view.  (In  our  simulation,  satellites  were  selected  from  prior  know¬ 
ledge  of  the  flight  profile).  This  would  also  require  building  into 

the  computer  simulation  programs  the  capability  of  multiple  reset  seg¬ 
ments  (such  as  in  SA?TS) . 

5.  Another  area  of  Interest  would  be  comparison  of  sequential 
versus  simultaneous  measurements.  For  example,  a  sequential  measure¬ 
ment  scheme  might  take  single  measurements  from  satellite  number  one 
at  t  •  15  seconds,  satellite  number  two  at  t  -  30  seconds,  satellite 
number  three  at  t  •  45  seconds,  and  satellite  number  four  at  t  ■  60 
seconds.  This  sequence  would  then  be  repeated  throughout  the  reset 
segment.  The  performance  of  such  a  filter  could  then  be  compared  with 
one  such  as  in  this  thesis  where  measurements  and  updates  are  performed 
with  four  satellites  simultaneously.  This  comparison  is  of  interest 
primarily  because  of  the  increased  cost  of  equipment  required  to  per¬ 
form  simultaneous  measurements  and  computations  as  compared  to  the 
sequential  case.  . 

6.  Various  flight  profiles  (for  example  a  "figure  8")  should  be 
simulated  to  Investigate  the  dynamic  effects  of  climbing,  driving  and 
turning  upon  the  filter  design  of  this  thesis. 
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Appendix  A 

Ten  State  Filter  Plots 

Plots  of  the  PXS  values  of  the  nine  IKS  plant  error  states  are 
presented  in  Figures  51  through  59  of  this  appendix,  The  structure  of 
this  10  state  filter  is  detailed  in  Chapter  VI,  These  plots  illustrate 
the  effect  of  elininating  too  many  states  from  the  filter  state  variable 
vector 
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Appendix  B 

Plots  of  the  Plant  Error  States  at  Specified  Update  P.ates 

The  plots  contained  In  this  appendix  are  the  remaining  six  plant 
error  states  at  the  various  update  rates  from  Chapter  VII.  The  plots 
on  the  following  pages  are  as  follows:  y  position  error  with  update 
rates  of  5  seconds,  15  seconds,  60  seconds,  and  90  seconds  are  shown, 
in  Figures  60  through  63,  Altitude  error  with  the  same  rates  are 
shown  in  Figures  64  through  67,  x  and  y  velocity  error  in  Figures  68 
through  75,  x  attitude  error  in  Figures  76  through  79,  and  azimuth 
error  in  Figures  SO  through  83. 
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